STUDY OF RANDOM FIELD MODELS 

FOR TEXTURE 


A Thesis Submrtted 
in Partial Fulfilment of the Requirements 
for the Degree of 

MASTER OF TECHNOLOGY 



By 

POOLAKMOY BANDYOPADHYAY 


to the 

DEPARTMENT OF ELECTRICAL ENGINEERING 

INDIAN INSTITUTE OF TECHNOLOGY, KANPUR 


JULY, 1983 



1 3 


» li Cd f 

twik'S ^ w" ■ itfj j«C ;I|i^^■- ■* 


6£- HSl- M-&AN-STU 




Certified that this work * STUDY OF RANDOM FIELD 
MODELS FOR TEXTURE* by POOLAKMOY BANDYOPADHYAY is 
carried out under my supervision and is not submitted 
elsewhere for a degree. 



{ S.K. Mullick ) 
Professor 


Department of Electrical Engineering 
Indian Institute of Technology 
K/iNPUR 


POST CH ' ‘ 
This ■ 

(OT ii tC HW3, ru : 

I 

I ill ‘ ‘ ' 

\ y . 

i,*..- ■ ’'Vr 


;‘.p* p)ved 
vec ut 
-.1.' Tech. I 


i i 3 T 1 



ACKInIOViJLEDGEMENT 


I vvish to take this opportunity to express my deepest 
sense of gratitude to my supervisor Dr, S.K. Mullick who 
has initiated me into this problem and provided me with 
infall iable guidance. His sincere advice, keen interest 
and constructive criticism during the course of this work 
has been a source of constant encouragement to me, 

I don't have words to express my gratitude to Prof, R,L, 
Kashyap. Long discussions with him during his stay in IIT, 
helped me a lot to take up the challenges in the field of 
image processing. 

I am also indebted to Prof. V.P, Sinha and Prof, J,N, 
Kapur for their constant encouragements and suggestions. 

Thanks are due to Mr, K.N.H, Bhat, P.G, Poonacha and M«K 
Patra for the fruitful discussions in various stages of 
the thesis, 

I am also indebted to DRDO for financial support. I also 
thank all others who made my stay at IIT memorable and 
pleasant. 

Finally I deeply appreciate the excellent typing done, 
by Mr, J,S. Rawat, 


P. Bandyopadhyay 



LIST OF CONTENTS 


Page 

CHAPTER 1 INTRODUCTION 1 

1*1 Texture Models 3 

1*2 The problem and the approach 4 

1.3 Literature Survey 6 

1.4 Outline of the Thesis 11 

CH/J>TER 2 CAUSAL SPATIAL MODEL FOR TEXTURES 13 

2*1 Introduction 13 

2.2 Review of ID time Series Model 14 

2.3 Principle of Parameter Estimation 17 

2*4 2D Statistical Model 23 

2.5 Empirical Equations for Model 

Building 30 

2.6 Model Building 32 

2.7 Simulation Results 35 

CHAPTER 3 FREQUENCY DOJ^^/illNl CH/\R/.CTERIZATION OF 

STOCH/iSTIC TEXTURES 42 

3.1 Introduction 42 

3.2 PN Sequences and /grrays 43 

3*3 Pseudo-random Arrays: Generation 43 

3*4 Spatial Filter Model 51 

3.5 Texture /uialysis 52 

3.6 Parameter Extraction 53 

3.7 Use of Separable Filter to Compute 

Phase 60 

3.8 Simulation Results 62 



Page 

CH/^TER 4 ANALYSIS AND SYNTI-ffiSIS OF TEXTURES BY 

NONCAUSAL RF MODELS 67 

4.1 Introduction 67 

4.2 Notations 68 

4.3 Summary of Discrete Spatial Interaction 

Models 70 

4.4 Infinite Lattice Models 72 

4.5 Finite Lattice Models 75 

4.6 Analysis of Finite Lattice Models 79 

4.7 Properties of FL Models 85 

4.8 Relationship between SAR and CM Models 87 

4.9 Parameter Estimation 89 

4.10 Choice of Neighbourhood for SAR and 

CM Models 93 

4.11 Experimental Results 97 

4.12 Discussion 106 

CHAPTER 5 TEXTURE SEGMENTATION BY SAR MODELS 107 

5.1 Principle of Segmentation 107 

5.2 Application of RF Model for Segmen- 
tation 108 

5.3 Principle of Clustering 110 

5.4 Simulation Results 115 

CH/\PTER 6 CONCLUSION 124 

REFERENCES 126 



ABSTRACT 


This study considers statistical modelling of a special 
type of images, called textures. Both causal and noncausal 
neighbourhoods are considered. We consider the 2D time 
series models for textures and the parameter estimation 
schemes for such models* We also discuss the frequency 
domain characterization of textures and application of PN 
arrays for this simulations. Among the non-causal models, 
we present two important ones; the simultaneous auto~ 
regressive (SAR) and conditional Markov (CM), Each of these 
models is characterized by a set of neighbours, a set of 
coefficients and a noise sequence of specified properties. 

To avoid excessive computation, we use an approximate maximum 
likelihood scheme for estimating the parameters. Asympto- 
tically consistent decision rules are given for choosing 
an appro piate SAR or CM model. Model studies have been 
substantiated using synthetic textures. An image segmen- 
tation scheme based on SAR model is proposed and studied 
in detail. The programs generated for this study are also 
useful for images other than textures. 
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INTRODUCTION 

The subject of image modelling involves the constru- 
ction of models or procedures for the specification of 
images* These models serve a dual purpose that they can 
describe images that are observed and also can serve to 
generate synthetic images from the model parameters. We 
will be concerned with the models of a specific type of 
images, the class of models for representing textures in 
image. Proper modelling of image has become important 
because, the efficient design of computer algorithms for 
image analysis, processing and synthesis can only be done 
using the framev^ork of an image model. 

There are four important areas of image processing in 
which texture plays an important role: classification [l» 23 » 
image segmentation [3], realism in computer graphics [4] 
and image encoding [5], Besides being an intrinsic feature 
of realistic objects, texture also gives important information 
on the depth and orientation of an object [6], Julesz [7] 
considers the problem of generation of familiar textures, 
an important one from both practical and theoretical view- 
points, Understanding texture is also an essential part 
of understanding of human vision [S], These considerations 
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have led to an increased activity in the area of texture 
analysis and synthesis. 

There is no universally accepted definition for 
textures. Part of the difficulty in giving a definition of 
texture is the extremely large number of attributes of 
texture that we would like to include in the definition, Vte 
consider a texture to be stochastic, possibly periodic, 2D 
image field. The most important attributes of texture 
are coarseness, contrast,^ directionality, line linkness, 
regularity and roughness. 

Most texture researches concerning texture modelling 
are characterized by the underlying assumptions made about 
the texture formation process. There are two major assump- 
tions, and the choice of these assumptions depends primarily 
on the type of texture to be considered in the study. The 
first assumption which is called the placement rule viewpoint, 
considers a texture to be composed of primitives. These 
primitives may be of random or of deterministic shape, such 
circles, polygons or even dot patterns. Macrotextures 
have large primitives, vmereas micro textures are composed of 
small primitives. The texture image is formed from the 
primitives by placement rules which specify hov/ the primi- 
tives are oriented, both on the image field and Vvith respect 
to each other. Examples of such textures include tiling 
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of the plane, cellular structures and a picture of brick 
wall* The second viev/point regarding texture generation 
processes involves the stochastic assumption. The place- 
ment rule principle for textures may include a random 
aspect in the stochastic point of view, however, we may 
consider that the texture is a sample from a probability 
distribution on the image space. The image space is 
usually a square lattice and the value at each grid point 
in the lattice is a random variable (non— negative) , Textures 
such as sand, grass and water are not appropiately described 
by a placement rule model. The key feature of these images 
is that the primitives are random in shape and cannot be 
easily described, 

1.1 TEXTURE NriDELS 

By the model of a texture, v/e mean a mathematical 
process which creates or describes the textured image. The 
primary object of texture modelling is the descriptixjn 
of real textures, A secondary goal of texture modelling is 
classification of textures. The numerical parameters of the 
model can be used as features to classify the textures. 
Texture modelling is an important concept in the sense that 
it gives the inherent features of the image and hence helps 
further processing, A survey of the well known and commonly 
used textures model schemes has been given by Ahuja [9], The 
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problem of characterization of stochastic texture is basically 
a system identification problem v/here the system is a two 
dimensional one, 

1,2 THE PRO ELBA AND THE APPROACH 

In this thesis, we consider different approaches for 
characterizing a stochastic texture. Among the various 
methods of modelling stochastic texture, we only consider 
the random field (RF) models as they are mathematically 
tractable. In RF nrodel of texture, we assume that the texture 
is generated by a linear system [lO] v/hen the system is 
excited by white noise. The different pixel values of the 
texture follow a probability distribution depending on the 
system and the probability distribution of the input. This 
problem of system identification can be tackled in two 
ways 

i) Spatial Domain Approach 
ii) Transform Domain Approach 

In spatial domain approach, we represent the input-output 
relationship by a difference equation. Based on the depen- 
dency of any pixel on other pixels in a given region, we 
define a spatial model to be either causal or noncausal 
and they are shown in Fig, 1,1, 
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Causal Model Non-causal Model 

X at point of interest 

Fig, 1,1: Spatial Neighbourhood for images 

In case of noncausal model, v;e assume that the value 
of a pixel depends on all the neighbouring pixel values in a 
given region. In spatial domain characterization, the 
problem is to estimate the coefficient or the parameters 
of the governing difference equation. Different statistical 
schemes are generally used to estimate the parameters. For 
noncausal model, we assume the image to be torus, finite 
lattice image and depending on the input sequence, these 
noncausal models are of two type 

i) Simultaneous Auto Regressive (SAR) 
ii) Conditional Markov (Q^) 

In SAR model, we assume that the input to the system is 
uncorrelated but it is correlated one in case of CM model. 
These types of noncausal models are very useful in texture 
analysis due to their realistic assumption about the 
neighbourhood. 
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In the transform domain characterization, texture is 
vievi/ed as the output of a homogeneous spatial filter excited 
by white, not necessarily Gaussian noise. Different stati- 
stical measures give the amplitude and phase response of the 
filter. In identifying the system, this method is not as 
good as the parameter estimation scheme but this is a 
computationally efficient schen® due to fast transform 
algorithms, 

1,3 LITERATURE SURVEY 

The random field model of textures has been studied 
by many researchers in the recent years. Depending on the 
type of image, the random field model has two subclasses 

i) Global models 
ii) Local models 

Global models troat an entire image as the realization of 
a random field. Hunt [11,12] has proposed a nonstationary 
gauss ian RF model for textures. He has shown the appro- 
piateness of this model by subtracting the local ensemble 
average from each image point and showing that the resulting 
image fits a stationary Gaussian model. In a recent paper, 
Hunt [13], discusses three kinds of nonstationarities; 

Case 1; Hons ta tip nary mean, nonstationary autocorrelation. 
Case 2: Honstationaory mean, stationary autocorrelation^ 
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Case 3: Stationary mean, nonstationary autocorrelation. 

The nonstationarity in the image mean is described 
by the array of individual means taken over a specified 
neighbourhood about each image point. For the autocorrelation 
function, the breakdovm of stationarity is related to the 
way the correlation function changes over the image. Hunt 
uses three attributes of the correlation function 

to describe its spatial dependence: 
its energy, its width, and its shape, Trussel and Kruger 
[14] claim that the Laplacian density function is better 
suited for modelling high pass filtered imagery than the 
Gaussian function. They have also shown that the basic 
assumption which allows the Gaussian model to be used for 
image restoration purposes are still valid under a 
Laplacian model, 

Nahi and Jananshahi [15] suggest modelling the texture 
by background and foreground statistical processes. They 
assume that the two portions are statistically independent 
random processes with known (or estimated) first two 
moments . 

Recursive solutions based on differential equations 
are coimnon in one-dimensional signal processing and have 
been generalized to two dimensions, Jain [16] investigates 
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the applicability of three kinds of random fields to the 
image modelling problem, each characterized by a different 
class of partial differential equations (PDEs): hyperbolic, 
parabolic and elliptic, A digital shape is defined by 
a finite difference approximation of a PDE, The class of 
hyperbolic PDE are shown to provide more general causal 
models than autoregressive moving average models. For a 
given spectral density function, parabolic PDEs can provide 
causal, serai-causal and noncausal representations. Elliptic 
PDEs provide noncausal models that represent two dimensional 
discrete Markov field and can be used to model both 
isotropic and anisotropic imagery. Jain argues that PDE 
models is based on a well-established mathematical theory 
and, further, that there exists a considerable body of 
computer software for numerical solutions. System iden- 
tification techniques may be used for choosing the PDE 
model for a given class of image, 

Pratt and Faugeras [17] and Gagalov<?icz [18] view 
texture as the output of a homogeneous spatial filter, 
excited by v^hite, not necessarily Gaussian noise, A 
texture is characterised by its mean, the histogram of 
white noise and the transfer function of the filter. For 
a given texture, the model parameters are obtained as 
follows 
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i) The mean is readily estimated from the texture, 

ii) The autocorrelation function is computed to 

determine the amplitude response of the transfer 
function. 

iii) Higher order moments are computed to determine the 
phase response of the filter. 

Inverse filtering yields the white noise image and 
hence its histogram and probability distribution. The 
inverse filtering or decorrelation may be done by simple 
operators. For example, for a first order Markov field 
decorrelation may be done by Laplacian operator [19], The 
whitened field estimate of the independent identically 
distributed noise process v/ill identify the spatial process 
in terms of autocorrelation function. 

If neither knov/ledge nor hypotheses about the type 
of random field are available, a class of local models is 
used. These models assume relationships among gray levels 
of pixels in small neighbourhoods. The modelling process 
consists of choosing the relationship and evaluating its 
parameters. Two basic categories of local inode Is arise 
from the joint and conditional probability formulations 
of variations in a neighbourhood proposed by Whittle [20] 
and Bartlett [2l], Besag [22], however described a 
number of difficulties with these approaches. He has shovm 
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that the constraints on the conditional probability 
structure are so severe that they actually control the 
particular nodels* Nevertheless, the conditional proba- 
bility approach serves as the basis for a eommonly used 
class of models, called Markov image models. Different 
choices for the set of a pixel’s neighbours give rise to 
different Markov models* The specification of the proba- 
bility distribution of a pixel’s gray level, given the 
gray level of its neighbours, defines a strict sense 
Markov field representation. Two of the more interesting 
models of this class have been discussed by Besag [23} and 
they are called autobinary and autonormal schemes. 

A close relative of the autoraodel is the simult- 
aneous autoregressive model( S/iR)*Ka.y5hap [24] notes that 
the symmetry requirements of the SAR. are automatically 
fulfilled without the need to place prior restrictions 
on them, 

Kashyap [25] uses circulant (block) matrices to 
define SAR models for an infinite array (array obtained 
by repeating a given finite image). He has also given 
the scheme for estimating the parameters of the model 
and decision rules for the choice of optimum neighbours, 
Kashyap also discusses all. these schemes for a doubly 
periodic array defined over a finite lattice. 
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The local models discussed so fax relate the gray 
level of a pixel to the gray level of its neighbours, 
using a conditional probability formulation. Through this 
approach has got considerable attention, its primary 
difficulty, hov^ever, is the high dimensionality of the 
joint probability densities, even for small neighbourhoods, 
making parameter estiiriation complex and cumbersome, 

Rosenfeld and Troy [26] and Haralick [27] suggest 
the use of two dimensional spatial dependence of gray levels 
for fixed distances and angular perations, Haralick and 
numerous other investigators apply features derived from 
the co-occurence matrix to various texture classification 
and discrimination problems. 

The entire area of syntactic image model is not 
covered here as it is outside the scope of this thesis, 

1.4 OUTLINE OF THE THESIS 

After a brief introduction in Chapter 1, in Chapter 2 
we have studied the 2D time series model of textures. 
Synthesis and parameter estimation schemes for different 
auto-regressive and moving averages processes have been 
discussed, A procedure for checking the adequacy of the 
estimated parameter is also given. 

In chapter 3, we have studied the transform domain 
characterization of texture. Instead of computer generated 
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random number, simulation is carried out by using non-binary 
pseudo-random array, being ideally white in the cyclic sense. 
The amplitude and phase response have been estimated for a 
separable spatial filter. 

Chapter 4 discusses noncausal nx>dels for a finite 
lattice image. These images are assumed to be toroidal* 

Both simultaneous autoregressive (SAR) and conditional 
Markov (CM) images have been studied. For parameter esti- 
mation, we have followed Kashyap* s method for obtaining an 
approximate maximum likelihood estimate and a decision rule 
for optimum neighbourhood has also been studied* We have 
also discussed the spectral density characterization for 
such models* 

In Chapter 5, we have proposed an application of 
spatial interaction model to texture segmentation. Simula- 
tions have been carried out for different segmenting curves 
and for all such cases the proposed scheme worked satis- 
factorily. Some of the disadvantages and error performance 
of this method have also been discussed. Segmentation 
scheme is studied only for SAR models. 

Chapter 6 sums up the contribution of the thesis with 
recommendation for further work along the same line. 



CHAPTER 2 


CAUSAL SPATIAL MODEL FDR TEXTURES 

2.1 INTRODUCTION 

Time series modelling of any observed set of data 
is an old concept and is used extensively in statistics as 
well as engineering. These models are also being used 
successfully in social sciences and economics. Their success- 
ful application to ID data has led researchers to apbly 
such models to multidimensional dataj in particular to 2D 
spatial processes. In the area of image processing such 
models are currently being used in image modelling and sub- 
sequently, using these models, to tackle problems of image 
restoration, texture classification, image data compression, 
spectral analysis etc, [ 9 ]. 

In recent years, however, considerable efforts have 
been made to model a 2D data array (e.g, images) by models 
which will take care of the dependency of a point on its 
spatial neighbourhood rather than on its *past* as the 
concept of ’past and future* are not defined naturally in 
such cases. 

In this chapter we will study two dimensional time 
series models and their parameters estimation, but this 
2D analysis will be done only after a brief introduction to 
its one-dimensional counterpart. 
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2.2 REVIEW OF ID TIME SERIES MODEL 

Let a time series be represented by Y-i •• fY ■ * • • • 
and let ^03^^ represent an independent identically distri- 
buted (II D) sequence of zero mean and variance a ^ , Then 

(i) 

the time series may be modelled by 


y . = ^ + 0) . + • • 


CO 


= + S a. w. , 

k=o ^ 


( 2 . 1 ) 


with a =1 
o 


or y^ may be related to by the eqn, 


Yi = PiYi^i + ^2^i-2 ••• ^i 


= Z 
j=l 




+ 


( 2 . 2 ) 


Here the eqn. (2.1) represents a moving average (MA) process 
of infinite order whereas eqn. (2.2) represents an auto- 
regressive process of infinite order. From eqn, (2.1) we can 
always get a moving average process of finite order (say q) 
as 



q 

Z 

k=o 


akO) 


i-k 


+ \x 


(2.3) 
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Here |j, represnts the mean of the output process. 

Again from eqn, (2»2) we can get an auto-regressive 
process of finite order (say p) as 


1 

'■j = £ P-iYi 

1 3 i-J 1 


(2.4) 


Eqn. (2.3) and eqn. (2.4) will be denoted as I\AA(q) and AR(p) 
processes respectively. By combining these two, we get an 
auto-regressive-moving average process where 


y. = |i + S p .y* . + Z a, CO- 

X . - J X*" J - K X-** 




k=o 


k i-k 


(2.5) 


Eqn. (2.5) is often written as ARMA(p,q) process. 

Assuming a zero mean process, eqn. (2*5) gives 


y-+a,y. -, + ...+ ay. v. = b x.+b,x. . 

1^1-1 P i-P oil 1-1 


+ ... + 


( 2 . 6 ) 


Relationship (2.6) can be represented by an input-output 
LTI model given by Fig. 2.1. 






2D LTI- 


System 




Fig. 2.1; A time series model of Texture 
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Taking Z~transform on both sides of eqn. (2.6) 

"1 ' * *5 

+ • • •+ )Y( z)=( b^+b^z 4* • »<+ b 


where the Z- transform is defined as 

oo 


Y(z) = Z 

is—oo 


-1 


So the transfer function of the LTI system is 


^ -k 
b + S b, z 


1+ Z 3i,z 

k=l ^ 

for one sided Z- transform. 


-k 


Case 1; 


Let us assume bj^ = 0 for k = 


1,2, . . 


yields 


G(z) = 


b 


1 + Z a. z 
k=l ^ 


By taking inverse Z-transform 


( 1 + Z a. z"'^ ) Y(z) = b.X(z) 
k=l ^ ° 


^i 


Z ( aic^'^i-k’^^o^i 
k— 1 


V 


we get^ 

qZ’'‘)X(z) 

(2.7) 

( 2 . 8 ) 

Eqn, (2.8) 

(2.9) 

(2.10) 


and 


( 2 . 11 ) 
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which is an auto-rogressive process of order p. So it is 
seen that in case of an auto-regressive process, the system 
is an ali-pole system. 

Case 2: 


If we assume = 0 for k = l,2,..^p, equation (2.8) 
reduces to 


G(z) = + 


q 

E 

k=l 



( 2 . 12 ) 


From eqn, (2.12) we get 


Y(z) 





q 

E 

k=l 



(2.13) 


and 

q 

Yi = S b X (2.14) 

This, in comparison to eqn. (2.3) is a moving average 
process of order q. 


Again it is seen from eqn, (2.13) that, in case of 
moving average process, the system is an all zero system. 
From eqn. (2.8) we see that the system in case of A3MA(p,q) 
process is a pole-zero system. In this respect ARMA comes 
out to be most general one. 

2.3 PRINCIPLE OF P,y^AMETERS ESTIMATION 

The auto— covariance for the case of AE(p) process 

is, 

Vie = 


(2.15) 



18 


where E[. ] represents statistical expectation. Then from 
eqn, (2,4) we get, by taking expectation on both sides 

Yjj. = ^p^k-p 

The autocorrelation coefficient is defined as 

Pk “ '^k^'^o 

Hence from (2.16) and (2*17) we get, 

k “ ^lf*k--l ^2fk-2 Pp'fk-p 

For k = 1,2,3, ,, ,p, from eqn. (2,18), we get a set of linear 
equations 

fl “ ^ 1*^2 fl ••• + Pp Pp-1 

P 2 " * ^2* ^pfp-2 

# 

p " ^if^ p-l'^^2'^p-2'^ ***'^ ^p 

This is the well known Yule-Walker 028 ] equations for 
AR processes. We can always find a matrix— vector relation- 
ship for eqn. (2,19), There is a fast algorithm developed 
by Levinson to solve such a set of linear equation for 
estimation the parameters [29], 

If we assume t)e the ith coefficient in a kth 

AR process, then Yule Walker equation becomes 



(2,16) 


(2,17) 
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^ ^2 **• fk-] 

1 ^ fi ••• fk-: 


fll 


( 2 . 20 ) 


fk-2 fk-: 


The left hand side matrix in eqn. (2.20) is of Toeplitz nature, 
so any fast method can also be used for solving such equa- 


tions,* In the set of coefficients ? 


i— 1 , * . , k 


, the last 


one, i,e,,0j^j^ is of great importance and it is called 
partial autocorrelation, and in an AR process gives an 

estimate of the order of the process. From eqn. (2.20) we 
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In fact, partial autocorrelation function of any 
process reveals an interesting property of the model* It 
can be shown that PAC (partial autocorrelation) function 
^kk AR(p) process follows 



for k 

< P 

‘^kk = ° 

for k 

> P 


Hence the PAG of a pth order AR process has a cutoff after 
lag p* 

The moving average process, in comparison with AR 
process# is more complicated to analyse. Let us take the 
MA(q) process of eqn* (2*3) for the sake of analysis 





( 2 . 21 ) 

Again, for such MA(q) process, the covariance is defined as 


Yk = E [yiYi.k] 

. ' . Vo = E [ypi ] 

= +... + ^ 



( 2 , 22 ) 

2 2 2 

where D = (l+a, + 0:2 + .*.+ a ) 

m 

2 

and cr = E[w.w.]= input variance 
tx) ^ X X 

(2.23) 
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Since we have assumed that the input sequence is zero mean 
and II D > it follows that 

E['co^I = 0 

and E[£i).w.] = 0 for i J 

From eqn. (2*21) it follows that autocovariance at lag k 
is given by 

Vk = (“k + ••• + VkV 

for k = 1,2, . . . ,q 

=0 for k > q (2.24) 

So combining eqn, (2*22) and eqn. (2, 24), we get, auto- 
correlation 



“g-k^o > 


for k— 1 , 2 , . . , q 


= 0 for k > q 


(2.25) 


From eqn, (2. 25), for k = 1,2,3,,,., we get a set of non- 
linear equations as follows 

fl°q = ^q-l“q 1 

fa^q ^ “q-2“q J (2.26) 

^ D = 2a„ 

iq q q j 
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These equations can be solved for getting the 
Gstimates of the parameters. It is evident that^ for a 
lower order MA pro-cessf the : solution may be obtainable but 
when order of MA process becomes large, the solution becomes 
less tractable. 


An autoregressive moving average (ARMA) process is 
much more complicated than both AR and MA process in, the 
sense that it contains two sets of parameters 
In ARMA process, the parameter estimation is again done after 
computing the correlation coefficients. In case of ARMA(Ptq) 
process, these will be q autocorrelations 
whose values depend directly on the choice of q moving 
average parameters a as well as p autoregressive parameters 
g.. Also the p values Provide the 
necessary starting values for the difference equation 
Cf(B) = 0 ¥k >, q+1, which entirely determines the auto- 
correlation at higher lags. If q-p <0, the whole auto- 
correlation function = 0,1,2,...., will consist 

of a mixture of damped exponentials and/or damped sine waves, 
whose nature is governed by the polynomial Cf(B) and the 
starting values. If hov\ever q-p >_ 0, there will be a 
q-p+1 initial values 5^q-p follow this 

pattern. These facts, are useful in identifying mixed 


processes 
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Since parameter estimation in ARMA processes is not a 
simple procedure. Box and Jonkins have given a 

chart (P--520) which will give the initial estimate of the 
parameters for a ARMA(l,l) process, 

2.4 2-D STATISTICAL MODEL [30] 

In the previous section we have studied the analysis 
of different time series models, viz,, AR, MA and ARMA, 

Those methods can be applied to 2D discrete images only 
when the image array is converted to a 1^ vector by some 
scheme (like lexicographic ordering). If z.. is representing 
a particular pixel value in the (i,j) position of the lattice 
then, it can be expressed again as a weighted sum of the 
data and white noise sequence (IXD) within a particular 
region or neighbourhood, A portion of such a lattice is 
shown in Fig, 2.2, and all, the grid points are equally 
spaced. 
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Hence in 2 dimens ion^^ the moving average process 
is given by 

— 

(2.27) 

If we define the horizontal shift operator by H, then. 


h'^z. . = z. . 


(2.28) 


Similarly a vertical shift operator is defined by V where 


z . . = z . 

1,1 i-m,j 


(2.29) 


Two more parameters need to be defined for the sak.e of 
analysis they are horizontal difference operator and 
vertical difference operator defined as follows. 

Horizontal difference Operator: 

Vk = 1 - H 


- V 


, 2. . = (l— H) z. *= Z. * — 2. . n 

h ^ 1,1 1,1 1,1 1 , 1-1 


and vertical difference operator: 


(2.30) 


Vv = 1 - V 
Vv ^1,J = 


i,l 


z. .-z. - . 

1,1 1 - 1,1 


(2.31) 


These operators are commutative 
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z. ^ 


m ^ n 


v; V 


V 


V h ^i,J 


If we define 


-p OO CJC» 

^ ^ n£o n"o V 


(2.32) 


then eqn. (2.27) reduces to 


Z.^J =^(H,V)U).^J With =;l (2.33) 

Since 5 '(h,v) represents the system whqse output is the 2^ 
time series|^z^ hence it may be called the system function 
of the model. 

Again in two- dime ns ion, the auto-covariance at a 
lag (kj£) is defined as 




(2.34) 


From eqn. (2.33) we get. 


Vk,p = E[‘i'(H,V)a..^ .$(H,V)a.j^^^^j^^] 


OO OO C30 oo 


z Z Z Z ^ti)^ 
p=o ■ •q=o m=o n=o P»‘I 


^^“i'-Prj-q ^i+k-ra,j+4^- n^ 
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This expression is nonzero only w^ien psin-k and cf=n-iL and 
for th^t value 


E[(jO, . f Oi. . ] = o 


oo oo 


= ^o4o'f 

The autocovariance coefficient y, p at lag (k,£) can also 

“T 

be derived from the model system function Let us 

define a function PCHjV) 


P(H,V) = r, S Yk 

' k=€6o 


( 2 . 36 ) 


Then from eqn, (2.35) we get, 

PCH.V) = <5^2 

= „ V“-P h"-'’ 

03 3 ^ ^ ni j n 

= ( SE V® h'^)(ES Vp H“^) 

w '■ Jm,n JP^q 

= ‘^(H,V) $(H~^, V"*^) (2.37) 

Again for the sake of analytical simplicity we study 
the estimation technique for 2D AR and MA process. 
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Case 1: 


AR(m,n) process: 

Here ^HjV) z. . = w. . 

if J if J 


where 




m n 


.ku£ 


k=o 


and d ^ = 1 
O f o 


Case 2; 


(2,38) 


MA(p,q) process: 


Here z^^^ = ©(H,V) 


P q _ £ 


where ©(H,V) = E I 9 ^ 

k=o 4=0 


H 


(2.39) 


Case 3 « 


ARMA (ni,m; p,q) process: 


Here ^(H,V) z^^^. := ©(H,V) 


(2.40) 


J.ike one dimensional AFIMA process, the two dimen- 
sional ARMA process is complicated to analyse and the 
estimation of parameters in this case is not a straight 
forward procedure. In the remainder of this chapter 

we will study the AR and MA process in two dimension and 

will give an empirical rule to estimate the parameters. Test 
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for stationarity is also given for a given discrete 
data set. 

Let us consider an AR( 1,1 ) and MA(l,l) model which 
are the simplest non-trivial 2-D processes. 



^1,0 ^^1,!^ ^0,1 ^k,^-l **■ ^l^l”''^k-l,^-l 

similarly 

;^k,^= ^l,of k-l,^L‘^*^0,lf k,t-l *^^l,lf k-l, f;-! 

Since we are considering AR(^,1) process, hence by putting 
pairwise value of (k,l) as (0,1), (1,0) and (.1,1) we get 
the corresponding 2D Yule-Walker equation as 
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The two dimensional correlation coefficient at lag (k,iD 
defined as 

?k,£ "" "^ 0,0 

The given matrix vector equation may be solved for 
the initial estimate of the parameter set. It is evident 
that dimension of the matrix and vectors depend on the order 
of process. 


Case 2 : 

MA(1,1) process* 


Here z. . = 0(H,V) w . . 

-L J XJ 

0(H,V) = 1 + + e^^HV 

Since the autocovariance of z.. at lag (k,i2) is given by 





from the above equation we get 

Yqo = (i + ®10 + ®01"^ ®11 ^ 

^10 = ^®10 ® 11 ® 01 ^ 

"'^01 ~ ^®01 ® 11 ® 10 ^ ‘^( 1 ) 

■^11 ^ %1 



So autocorrelation functiohs are 



^10 

«10 + ®01 + ® l ! 


fox = 

fxx = 


^01 ^ ^ 1^10 
+ ®01 + 


9 


UL 


2 2 2 
1+ 9,^ + 9qj_ + 9;^i 


UO 
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Frofii these equations we can solve for the parameters 9 ^q# 
9qj^ and 9^^. One thing to be noticed here is that the 
complexity of the parameter estimations increases with the 
order of the model, but in case of AR process if the order 
of the model is increased then the method for solving them 
remain unchanged in the sense that v;e have to solve more 
number of linear equations and same Yule Walker equation 
can be applied* So given a set of discrete data v^ether it 
is one .dimensional or two dimensional, it is a general trend 
in all fields (where time series model works) to search 
for a suitable AR model. 


2.5 DAPIRIGAL EQUATI.ONS FOR.iMDDEI, BUILDIN3_ 

In application of these methods, to a given image 
data, the ensemble averages will be replaced by sample 
averages for computational conveniences. 



31 


For an image 5 where its dimensions is x n 2 » 

we define: L '^-3 


1) Mean 


_ X fell ^2 

2 _ ^ 2 2 z • ♦ 

”l 2 i=l j=l 


2) Autocovariance 


n 0 

te-k ’^2"^ 

Vk,^= Tn^-k) (n^-^ 4 ^ 


where (k,t) is the lag in (i,j) direction. 


3) Variance 



4) Autocorrelation 

5) Partial Autocorrelation function 
a) In horizontal direction 
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b) In vertical direction 


fkk = 


0,1 


k = 1 


f: 


A -P 

k,o ~ Y(k-l),jT(k-j),0 

k-1 

1 - Z Ct5 

j=l 


k = 2,3,..* 


(k-l),j j,o 


where 




ct, 


Ck-2),j 'fdc-D.Ck-D^Ck-a), (k-j-1) 

■Yj = 1,2, ...(k-1) 


2.6 MODEL BUILDING 

Given a set of 2D data we can not say by looking at 
it what type of model will be best for it. The information 
about the probable model (whether AE, MA or ARMA) can be 
obtained from the partial autocorrelation function. Box and 
Jenkins have shown [28] that the partial autocorrelation 
function gives an idea about the nature of the model as well 
as its probable order. They have given a chart of the 
nature of partial autocorrelation which is given here for 
AR, MA and ARMA. 

We have also discussed in the previous section that 
a data set can be modelled as a time series nxjdel if and only 
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if it is a stationary one, A data set is said to be 
stationary if its covariance matrix is positive definite. 

The simple way to check stationarity of a series is by 
the behaviour of its autocorrelation function. In a 
stationary sequence the autocorrelation decays uniformly 
without any notch. In case of 2D data we have to test the 
smoothness of the autocorrelation in horizontal and vertical 
directions. With these things in mind we can give successive 
steps for model building for a given 2D data* 

Step 1 : Calculate the mean value 

Step 2 ; Calculate the auto covariances 

Step 3 ; Calculate the variance of the data sample 

Step 4 ; Calculate the autocorrelations 

Step 5 ; Check stationarity of the given data, A given data 
may not be stationary in vertical or horizontal directions. 
This behaviour will be obvious from the autocorrelations. 

In most of the cases, fortunately ^it has been found 
that the difforoncos in the adjacent , pixel values in 
the nonstationary process remove the trend and make the 
process stationary. So, in an attempt ( in general) to 
remove nonstationarity, we take differences in horizontal 
and vertical direction of the given 2D data till a stationary 
data series is apparent. 
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Take the difference in the direction in v;hich 
the data is nonstationary. The information about the 
direction of non-stationarity- can' be oMainecL from, th® 
p a r tia-i auto co r e la t io n f un c t io n * 

Step 7 ; Determine the partial ACF in order to check the 
approximate order of the model and check stationarity. 

Step 8 ; Identify the nK>del from Table 2*1, 

Step 9 : Estimate the parameters. After identifying a 
tentative nwdel we can use the sample autocorrelations to 
obtain a preliminary estimate of the parameters. These 
values may be used as an initial estimate for obtaining maxi- 
mum likelihood estimate. 

For computing initial estimate of parameters in case 
of AR process, we have to solve a set of linear equations 
given by Yule Walker equation. But in case of MA parameter 
estimation, we have to solve a set of nonlinear equations 
and in that case determination of preliminary estimates is 
less tractable. 

After identifying a tentative set of model parameters 
for a given pictorial data, we v;ill obtain the best-estimate 
of the parameters. This can be done by checking the rando- 
mness of the residual. The concept of residual in time 
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Series model is as follows* Suppose, for any pixicess 
(whether AR or MA),the actual value of parameter vector is 
9, Given a IID sequence (preferably Gaussian), we can 
always generate an image y by using the IID sequence and 
vector 0. Let us say © is the estimate of the parameters, 

Vi?ith the Same HD sequence and parameter 4 if we generate 

. A A 

— an— image ^ then the difference between ^ and ^ is call'! 

called the residual. If the model is best one, then the 
residual = X “ i) will be white noise sequence. This 

is what will be achieved in step 10. 

Step 10 ; Compute the residuals and compute its autocorre- 
lation functions. The autocorrelation of the residual will 
suggest the direction in which the model should be modified. 

Step 11 ; Make modifications in the model based on the 
obtained information. 

All these steps are also given as a flow chart in Fig, 2, 3, 
2.7 SIMULATION RESULTS: 

For simulation purposes, the textures were first 
generated using specific models and then the parameters 
were estimated from such synthetic textures. 

For this purpose, we have only considered lower order 
AR and ^'lA processes. Programs are also developed for 
generation of synthetic textures and subsequent parameter 




Fig» 2#2s Modol Building Flow-«chart 
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estimation and they are given in Appendix I, Provisions 
for handling higher order models is incorporated only. for 
AR processes because the nonlinear nature of the autocorre- 
lation equations in M/v process cannot be predicted from 
the knowledge of lower order MA models. 

In this study of 2D time series models, we have 
generated several synthetic images of size 64x64 and 32x32, 
Some of such synthetic textures are shown in Fig, 2*4, In 
these cases we have used only AR(l,l) and MA(1,1) processes. 
Gaussian textures were only considered and for synthesis we 
have used IMSL.(Intn’l Mathematical Subroutine Library) 
routine GGNOR for generating Gaussian IXD sequence. But for 
2D time series nK>dels we need 2D, IXD sequence or an array. 
Since the elements of the sequence are independent so any 
arbitrary arrangement of them does not change their inde- 
pendence, Thus one way of generating a MxM uD array is to 
generate a M xl LID vector and then arrange them (lexicogra- 
phically or by any other means) to form a MxM array. But 
considering the space requirements, v/e took another approach. 
In this method, each row (or column) of such IID array, was 
generated by using different seed value of GGI^R routine. 
From here onwards whenever such computer generated random 
number array with prescribed density is generated we have 
adopted the same scheme. Though we have mentioned in the 
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last section that for best estimate of parameters we have 

to take the maximum likelihood estimate, but it has been 

found that for image size of 64 x 64 the initial estimates 

are quite close to the actual parameter values which is 

explainable from statistical . point of view. Instead of 

computing the maximum likelihood estimate which is a slow 

process in computers we have checked the whiteness of the 

residuals. Here we will give in a tabular form the 

actual and estimated parameter set for AR(1,1) and MA(l,l) proces 

process for a 64x64 image. We have also given the plot 

of the autocorrelation function of the residuals for finite 

number of points for AR(l,l) and MA(1,1) processes, in 

Figs, 2*5 and 2,6, 
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Table 2.1: Summary of the Properties of AR» MA and ARMA Process [30] 
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Table 2.2 

IMAGE SIZE : 64 X 64 
MODEL : AR(l,l) 

Input : GAUSSIAN (0,1) 


Actual Parameter Estimated Autoco- Estimated para- 

Values rrelations meter Values 


= 0.40 
C^Ol = 0.40 
^^11 = 0 -^ 


00 


1.000 


01 

10 


0.8551 

0.9299 

0.8214 


= 0.21 QQ = 1.0000 


=-0.14 

= 0.12 


01 - 


0.7216 


10 

11 


0.8126 

0.6924 


= 0.28 
= - 0.12 
= 0.14 



= 0.3601 
= 0.4612 

^11 ^ 

= 0.1832 

^01 "" - 0-1216 
= 0.1522 

= 0.2612 

^ 0.1 ” 0 ^ 10^3 

= 0.1567 


= 0.7964 


Table 2,3 
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IMAGE SIZE : 64 x 64 

MODEL ; MA(1,1) 

INPUT : GAUSSIAN 

(0,1) 

Actual Para- 

Estimated Auto- 

Estimated Para- 

meter Values 

correlations 

meter Values 

®10 = 0*^ 

00 = ^-0^° 

©10 = 0.3521 

= 0.40 

^0 = 0.34698 

®01 = 

. 

0^ = 0.34698 


9 ^^ = 0.40 


© , = 0.3523 


= 0.25327 


SlO = 

00 ~ 1*0000 

©^0 = 0.2723 

601 = 

= 0.2634 

©.., = 0.1412 
. ul 


0j_ = 0.1654 



= 0.1448 

= 0.1623 

©10 =,0.12 

... = 1.0000 
uv 

^ U # 

©,.-^ 1 = 0 • 14 

UX 

= 0.1686 

XU 

^ = U^-I‘u23 

vJ X 


31 = 0.1176 


^11 = 

= 0.1613 

* *3 • 1721 
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CHAPTER 3 


FREQUENCY DOMAIN CHARACTERIZATION OF STOCHASTIC TEXTURES 
3,1 INTRODUCTION: 

In the last chapter, we have studied the spatial 
domain characterization of textures by standard two dimen- 
sional time series models. It was also stated that models 
based on either autoregressive or moving average processes 
were computationally efficient. In general, parameter 
estimation for ABMA process is a complicated procedure. 

In this chapter, we shall study another procedure for chara- 
cterizing the texture. This procedure represents the system 
in terms of amplitude and phase response. This' approach 
may be thought of as a frequency domain equivalent of time 
series model. 


White 

noise 


Homogeneous 

Spatial 

System 




Texture 


Fig, 3,1: A model for stochastic texture 


The Figure 3,1 ^ows the model of texture generation 
Like the two dimensional time series model, here also we 
assume that stochastic texture is generated by passing a 
white noise sequence through a linear homogeneous spatial 
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system. In our study, the 2D system is a 2D digital filter. 

In literature, the white noise sequence is generally taken 
as the computer generated random number sequence or IID 
sequence, having some specified probability distribution. 

But, in the frequency domain characterization of textures 
we use a 2D pseudo random. array as the input noise array, 
pseudo random sequences (generated from shift register 
sequences )have the property that they show nearly ideal 
autocorrelation function when defined in a cyclic fashion 
and they are also periodic. So, when such a doubly periodic 
array is passed through a linear system, the output also 
becomes periodic. So computation of amplitude and phase 
response can be done by FFT techniques. The useful concept 
for such system characterization is the moments of different 
orders. This method of texture modelling was tested here by 
generating a synthetic texture by such a process. For 
this simulation, nonbinary PN array was used as noise input. 

.3,2 PN SEQUENCES AND ARRAYS: 

In this section, we will study the relevent properties 
of pseudo random sequences and its generation. It is found 
[31] that sequence, v\hen generated as a maximal length 
sequence, shows some good properties. Such binary sequences, 
if they are generated by a m length shift register, will 
be periodic sequences of period 2^~1. If we compute the 



44 


cyclic autocorrelation of such an n ler^th sequence (n=:2®'-i), 
the the autocorrelation follows the given relation 

"^(0) = 1 (3.1) 

j^(i) = - ^ , for 1 < i < n-l 

where ■^(,) is the normalized autocorrelation function. Such 
a binary PN sequence approximates a white noise sequence, 
because y for a large length SGquence,-^(i) j'^i ^ 0 becomes 
negligibly small, which leads to a periodic delta function 
structure of autocorrelation. Such sequences have been used 
extensively in digital data scrambling, synchronisation etc., 
and a considerable body of literature exists [32,33,34 3. 


Recently, several applications [ ] have called for two 

dimensional arrays whose 2D autocorrelation function is such 
that'^(0,0) = 1 and^(i,j) is small for (i,j) 4 (0,0}. It is 
very easy to generate such arrays from a given PN sequences 
[31]. For generating such an array of dimension nj^xn 2 with 


p(0,0) = 1 

’ 0 ^ <3.2) 

0 < j < n2-l 
(i,j) 4 (0,0) 


we need a sequence of length n, where n = n 2 ^n 2 = 2'^-l, 
provided n^^and n 2 are relatively prime. For generating a 
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31x33 array (as in our case), we need a PN sequence of 
length 1023, Since the generated sequence should be maximal 
length sequence in a specific Galois field [ ], we have 

to take the required primitive polynomial in that field. 
Primitive polynomial in a shift register sequence dGte3>- 
mine's the governing difference equation. For the 

generation of a n-length shift register sequence (n=2’^-l), 
we need m initial values to start with. So, given the 
initial values and a primitive polynomial in the specified 
field, we can always generate a maximal length sequence. 

The most important property of such maximal length sequence 
(MLS) is its autocorrelation function. The autocorrelation 
function '^(i) of a real (or complex valued) sequence s^, 

^1' ^2’ *** ®n-l length n is defined in cyclic sense as 



s .s 
J 




i — 0 ^ ^ 


(3.3) 


where bar represents complex conjugate and j(7) means 
summation modulo-n. This autocorrelation function J^(i) 
is a periodic sequence ^i.e, (n+i) = ^(i). The auto- 

correlation function is shown . in Fig, 3.2. The auto- 
correlation function of a pseudo-random sequence of length 
n = 2’®-! is given by 
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Fig* 3,2: Cyclic ACF function for a PN ^sequence of length n 

w «* Ol 

f (0) = 1 

j^(i) = - i 1 < i < 

Instead of binary PN sequence , having entries from 
0 and 1 only, we can also have a nonbinary sequence, taken 
from an alphabet of q symbols, for any q which is a prime 
or power of a prime. 

Let q be any prime power and let GF(q) be the Galois 
field with q elements. We need a polynomial 9(x) of degree 
m with coefficients from GF(q) which i) is not the product of 
two such polynomials of lower degree, and ii) has a zero as 
a primitive element of GF{q”^), say ^ , V/e call such a 
polynomial primitive over GF(q), In case of GF(4) (=GF(2^)) 

I 2 2 

we have the elements < 0,1, ? where m is a root of 
the primitive polynomial of degree 2 in GF(2), If we take 

such polynomial as x +x+l then, since m is root of it, 

2 3 

u + w + 1 = 0 and m = 



1 
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Suppose, the primitive polynomial is 

g(x) = + h^x+h^ (3.3) 

with h^^GF(q), h^ 4 0. The shift register specified by the 
above polynomial is as follows, in the Figure 3.3, the 



Fig. 3,3: General scheme for PN sequence generation 

boxes contain elements of GF(q), say a^. The 

feedback path then forms as 

^i+m mf“l i+n>“l nv-S^i+m—S ** * l^i+1 o^i 

(3.4) 

Eqn, (3,4) gives the recurrence relationship describing the 
output sequence. This is an infinite sequence of period 
q’^-1. A segment of such sequence of length q^^-l is called 
the pseudo random sequence over GF(q). As an example, if / 
q = 4, m = 2 and g(x) = x +x+w, we get the PN sequence 
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Olltu^lOiowlwOw^w^ww^ (3.5) 

of length 4^-1 = 15, This sequence can be modified to 
a sequence 

011310221203323 (3.6) 

by replacing w and w by 2 and 3 respectively. 

To find the autocorrelation function of a nonbinary 
PN sequence, we have to get a complex valued sequence 
from the real one in order to use eqn, (3.3). The unique 
transformation (from real to complex) which gives the 
autocorrelation as in binary sequence is as follows. 

If q is a prime or prime power then the complex 
sequence a is obtained from a by replacing each Y£_GF(q) 
by exp [Y^-1 After getting such complex sequence 

if we use the eqn. (3.3) then we get, 

^(0) = ^ 1 < i i ^"-2 

3,3 PSEUDO-RANDOM ARRAYS :GENERATION! 

Though there ard umpteen number of ways of making 
a matrix out of a vector, only one method out of them is 
useful in making pseydo random arrays so that the array has 
properties similar to that of the generating sequence. 

In the context of pseudo random arrays it is needed that 



49 


this 2D array will have same autocorrelation structure as 
that of pseudo ..random sequence. 

Firsts we will consider the binary sequence and 

corresponding array, A binary (modulo-2) PN array is 

constructed by folding PN sequence. Let us take number of 
k-,k« k, 

the form n=2 -1 such that n^ = 2 -1 and = n/n^. 

Examples are 

n=15=2^-l with kj_=k2=2 n^=3, n^ = 5 

n=63=2^-l with kj_=3, 1^=2 = 7, n2=9 

Cyclic autocorrelation of such array (nj^xn 2 ) is given by 

-{^(o.o) = 1 

f ^ (0.0) 

The condition to be fulfilled is that n^^ and n 2 are relati- 
vely prime. Here we will explain the construction procedure 
of such arrays. 

The starting point is a PN sequence a 
obtained from a primitive polynomial of degree m=K^K 2 . We 
use a to fill up an n^^ x n 2 array by writing a down the 
main diagonal and continuing from the opposite side whenever 
an edge is reached. Figure 3,4 shows the construction of 
such array from a PN sequence when n=15,nj_p=3 and n2=5* Thus 
the PN sequence 
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00010 0 110101111 
produces the array 


0 

1 

1 

1 

1 

0 . 

0 

1 

1 

0 

0 

1 

0 

0 

: 1 


(3,7) 


(3.8) 


Start at 
north-west 
corner and 
move south- 
east 



Continue from 
opposite side 
when an edge 
reached 



The complete 
array 


^0 

^6 

^12 



^10 

^1 

a? 

^13 

^4 

1 

^11 



®14 


Fig, 3.4: Construction of a PN array 
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These arrays are doubly periodic with period 
(horizontally) and n 2 (vertically). For such periodic 
array we can also define 2D cyclic autocorrelation function 
as 




(3.9) 


and 0,1> , • • y (n^”!) 


where x(i,j) represents the PN array and ^ is addition 
modulo n. 


For such a PN array, it is found that 


R(0,0) = 1 
and, R(k, 


(3,10) 


So 5 for large values of n^ and n 2 j the autocorrelation 
will be an approximate two dimensional delta function which 
will give an approximate flat power spectral density. 

Same procedure holds for generating non-binary sequences. 


3.4 SP/vTI/X FILTER MODEL 

In this section, we will try to model a texture as a 
homogeneous 

realization of a ^ spatial process. Though it is 


w 


/ « 
■ 5*3 • 
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a standard procedure to create stochastic texture from 
computer generated random numbers, in this thesis we have 
studied textures generated from PN arrays. We make the 
assumption that the process is homogeneous and ergodic. 

Though homogeneity is not a valid assumption for real images, 
it may be a good one for texture type images. The spatial 
process or the texture is the output of a homogenous spatial 
filter,- Most of the time series nx^dels assume that the 
white noise sequence is a Gaussian one, but here we do not 
make any such assumption. In fact, most of the real life 
textures are non— Gaussian, Since iMdelling of a random 
image means finding the system which will generate the 
image when excited by white noise, we can represent the 
system in two ’Ways, i) in sample domain and ii) in frequency 
domain. This particular model represents a texture (or 
homogenous image) in terms of the frequency domain character- 
ization of the system, 

3.5 TEXTURE ANALYSIS! 

In this section we will show that if a texture can be 
represented by the proposed model in Fig, 3,5, then it is 

^ Texture array 


Fig, 3,5: Proposed Model 


PN array 




Spatial Operator 
0 
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possible to compute all the filter parameters and then by 
inverse filtering we can get back the input noise 
parameters, 

a) Set of Texture Parameters: 

When a texture is generated or is supposed to be 
generated by the above nradel, then the resulting parameter 
sets for the textures are 

i) the texture mean 

ii) the histogram of the input white noise 

iii) the transfer function H(z) of the filter (or its 
impulse response). 

The second subset is equivalent to the set of different 

moments of the input IID sequence. These moments may bo 

defined as where 

i 2 n 

= nth order moment 

= E[b^(k)] 

where b{k) is the input noise. 

This set of moments and histogram of the noise are equivalent 
in that sense any one of them will completely characterize 
the noise sequence, 

3.6 parameter extraction [18] 

In this scheme of image mdelling, the process of 
parameter extraction can be given in block form as in Fig, 3, 6. 
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Fig. 3.6: Block diagram for parameter estimation 

3,6.1 Amplitude Response Estimation: 

Let us assume t(kj^,k 2 ) be the pixel intensity at 
the point (kj^,k 2 ). V/e also assume that input array 
b(kj^,k 2 ) is doubly periodic with period N1 and N2 in 
horizontal and vertical directions respectively. Let us 
assume that h(m,n) represents the impulse (2D) response of 
the filter. So given h(m,n) and b(k^iC 2 )»wo can write the 
output t(k 2 _#k 2 )as the 2D convolution of h(m,n) and 
b(ki,k2)as 

t(kj^,k 2 )= h(m,n) * * b(kj^,k 2 ) 
v^/here * ^ represents the 2D convolution 
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or, 

Nl~l N2-1 

tCk^k^) = b(m^,m2) h(k3_~m^,k2-ni2) (3.11) 

Here we can also define 2-dimensional Z- transform of an 
array x(i,j) as 

oo oo * • 

x(z, ,z.) = E E x(i,j)z z (3.12) 

i= -<o J=r — oo ^ 

So if we define T(zj^,Z2)» H(z^,Z2) and B(Zj^,Z 2) as the 2D 
Z- transforms of the array tCk^^^jk^), h(kj^,k2) and b(kj^,k2) 
then from eqn. (3,11) we get> 

T(Zi>Z2^ = H(Zi,Z2).B(zi>Z2) (3,13) 

Though the first order moment or mean of texture t{k^,k2) 
is defined statistically as 

t(kj^yk2) — 
it is estimated as 

_ 1 N^-1 N2-I 

^a(^1>^o) ~ M M ' ^ ^ (3.14) 

e 1 ^ Nj^N 2 r z 

The second order moment of the texture t(kj^,k2) gives 
the autocorrelation function array ^2^ texture. 

The second order moment of t(kj^,k2) is defined again in 
statistical manner as 
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M, 


= E[t(k^,k2)t(k^+!?^, k2+ 4>] (3,15) 

Using eqn. (3.11) we get 


M. 


Ni-1 

1 

CM 

E[ s” 

z 

0 

II 

m2==o 


Ni-i 

N 2 -I 

B E 

£ 

mj^=o 





(3.16) 




(3.17) 


Eqn, (3,17) is true only if input array is a white one. 

Here B 2 represents the 2nd order moment of the input 
noise. For time being, we can assume B 2 = 1* 

If we take the 2D Z- transform of we get 

00 M 

=:«-(» =2-^ ^ a. j. z. 

(3.18) 

H(Zi,Z 2 ).H(zi“^, 22 ^) (by using eqn, (3. 17)) 

(3.19) 

|H(z^,Z2)}^ (3.20) 

So we see, in 2D case, the 2D Z- transform of the 
second order moment of the output of a system gives the 
amplitude square response of the system. 
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If the process is Gaussian then the 2nd order moment 
and first order moment contain all the statistical infor- 
mation of the process. By inverse filtering with a zero— phase 
filter we can get back the white noise. In that case, it 
is not necessary to compute the higher order moments. But 
in case the process is not Gaussian, it is possible to 

estimate the phase of the filter from higher order 

moments, 

3,6,2 Phase Computation: 

Phase information is contained in the non-Gaussian 
part of the higher order moments. In the simplest non- 
Gaussian case we can assume that the process is not Gaussian 
at the third order i.e, B^ 0, We can define third order 
moment as 

M3(ki,k2jt.i,l2)=^C't^0,0)t(ki,k2)t(?i,22)] (3.21) 

It shows that the general filter structure gives a four 
dimensional third order moment function. In order to avoid 
long computation we describe the process of phase estimation 
by considering a special filter - called separable filter. 

If a 2D filter is a separable one, then it can be thought of 
as a combination of two one-dimensional filters. In case 
of separable filter, the 2D impulse response h(m,n) can 
be written as 

h(m,n) = hj^(m)h 2 (n) 


(3.22) 



58 


or in frequency domain, 

HCz^jZ^) = H^(z^)H2(z2 ) (3.23) 

where h^(m)<^:r^ H^( z^) 
tnd h^Cn) ^>H2 (z2) 

ViJe will show here hovj to use the concept of third order 
moment in. the estimation of phase of a unidimensional filter. 
For this purpose let us consider the Fig* 3.7. 

t(n) 


Fig. 3.7: One dimensional model 

Let us define the third order moment of t(n) to be 
M 3 (?^,? 2 ) where 

= E[t(0)tcii)t(£2) (3.24) 

t 

Mo(t,»2o)= Bq h(m)h(nH-{,)h(nH-?2) (3.25) 

be 

Here we assume b(n) and hence t(n),to £N~periodic. Taking 
Z“ transform of 

= 2 2 ^ ^2 ^ (3.26) 

^2 
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Using eqn. (3.25), we get 

= B3 H(z^)H( 22 )H [(z^jZ^)-^] (3,27) 

Let us assume z^ = e ■'■ and z^ - Q , then from 
eqn, (3,27) we get, 


+ Cf(w2)-C((wj,+a)2) (rnod 2ii) 


(3.28) 


where 


A 


^3^ ^1*^2^ == phase of M3(zj^,Z2) 


0( w) = phase of H(zj^) 

^ = 0 or 1 depending on sign of B2'. 

A 

So, after computing the phase of ^(^^(Zj^jZ^)# we have 
to find the phase of the filter H(z) i,e, 0(w) by using 
eqn, (3., 28). Here we will give an iterative procedure to 
compute the phase. 

Let us assume that we want phase 0(a)) of an array of 
N values, depending on the size of the filter. 


Let us define = 


2-n:(k-l) 

N 


and 



Hence from eqn, (3,28) we get 
03(k,?)= ii'H0 (k)+ 0(-£) - 0(k+4.-l) (mod 2%) 


(3.29) 
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The eqn* (3,29) is the digitized version of eqn. (3,28), 
For^= 2 in eqn. (3.29), we get, 

03 (k, 2 ) = 11 + 0(k)+0(2)-0(k+l) (mod 2 ii) (3.30) 

03 ( 1 , 2 ) = It + 0(l)+0(2)-0(2) (3.31) 


0(1) = 03(1,2) 


Similarly, 


assuming 1 


0 ( 2 ) 




(3.31) 


(3.33) 


We can rewrite eqn, (3,30) (with l) aS 

0(k+l) = 0(k) + 0(2) - 03(k,2) (3.34) 

Vi/e will use eqn, (3,29) to eqn. (3.34) for the phase 
estimation of the 2D filter. 


3.7 USE OF SEPARABLE FILTER TO CXDMPUTE PHASE 

Vile have already stated that if the linear filter is 
not a separable one then third order moment yields a four - 
dimensional function and from such function phase compu- 
tation is a very time consuming jprocess. For this reason, 
we have taken a separable filter. In such cases, the 
convolution also becomes very convenient one. 
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A filter h(i,j) is said to be a separable one if 
we can write 

h(i,j) = hji)h^(j) (3.35) 

in Z-domain H(z^,Z 2 ) “ (3.36) 

If we say X(z^,Z 2 ) is the Z~ transform of the input and 
Y(Zi,Z 2 ) is the Z-transform of the output, then entire 
process can be shown as in Fig. 3.8. 



Fig, 3,8: A Separable filter structure 

Though we can get Y by 2D convolution of X with H , 
still we can utilise the separable nature of. the filter. 
In the Fig, 3.8 the filters Hq and are both ID filter. 
Let us assume X is a n^ x array. Then the first filter 
Hq will take each column of X in succession and will modify 
them to get another intermediate array Z. . Each column 
of Z is the column filtered' version of X. This new array 

^ is fed to the next ID filter will take each row 

of Z and it will modify the row to produce a row of the 

array Y, So Y is the row-filtered version of the array Z,' 
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This separable property of the system can be exploited in 
computing the phase. If represents the phase of 

the ID filter and represents the ID phase of 

then phase of the total filter H is given by 

• (mod 2ii) (3.37) 


So effectively the 2D phase of the system is computed by 
the help of its one- dimensional portions, 

3.8 SIMULATION RESULTS: 


In the simulation process we have taken a 31x33 doubly 
periodic pseudo random array and a 2D filter with transfer 
function given by 


2-z, 2-Zo 

H(zi,Z 2) = (i.22 ) (i«2Zo) 


(3.38) 


This filter is a constant magnitude (all pass) filter. It 

has the property of introducing no correlations. By 

-jw, -jo) 

replacing z^ = e and ^2 = s we can compute the 

amplitude response of such filter on the unit circle. Then 




H(a)i,m2) 


l^ 




1^. 


-J(Or 






(l-2e •^) (l-2e ^) 


Since, e = cos - jsin co^ 

->2 • • 

3.n.cl ^ — cos ^2^ "" j s xn co^ 


(3.39) 


(3.40) 
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simple trigonometric identities lead to 
1 H f 10^ ) I = 1 - 


(3.41) 


This type of filter shows a pole- zero structure in its 
transfer function which in sample domain, is equivalent to an 
ARlvlA process. Since H(^) = v/here N(z) and D(z) are 

polynomials in z such filter is an HR (infinite impulse 
response) filter. But impulse response of such filter is 
also a decaying one in both direction. The input white 
noise array was taken to be 31x33 from GF(4), i.e. the 

elements of the array are from the set |^0,1,2,3^ . Since 

total number of elements in this array is 1023 (=4 -1) v/e 

need a shift register of length 5. V;ie also need a 

primitive polynomial of degree 5 in GF(4), The polynomial 
is taken to be g(x) 

g(x) = x"^ + X + 0 ) (3,42) 


where m is a root of the primitive polynomial of degree 2 in 

2 

GF(2). If we take that polynomial as x +x+l = 0 then 


2 

w + w + 1 = 0 
and = 1 


(3.43) 


So the required shift register is as shown in Fig. 3,9. 
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Fig. 3,9: The feedback connection for Eqn. (3,42), 

After generating the periodic sequence of length 1023 
we then used the scheme for converting a vector into a 
matrix to form a 31x33 array. This 2D array x(say) was then 
filtered by the spatial system H(z^,Z 2 ) to get the periodic 
output which Was the required texture. Though we have 
represented the filter in Z-domain by the transfer function 
Zz)* the filtering scheme was developed in sampling 
domain by writing the corresponding difference equation. 

We also took advantage from the separable property of the 
filter. Here we refer to the Fig, 3.8, For the given 
system 

2-z, ' 2-z^ 

^^^l’^2^ (^^22“) . (x.2z2 ^ 
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If we assume, 


H, 


= ^r=i::2z7 


(3.44) 


then, by taking Z- transform of eqn. (3,44), we get. 


y*(n)=G,5 x* (n)-x* (n-l)+0.5 y*(n-l) 


(3.45) 


As mentioned earlier, will take each column of the 
input and produces a new column following eqn, (3,45), If 
z represents the 2D output of then, 

z( i,3)=0,5x( i, j)~x( i, j-l)+0.5z( i, j-l) (3,46) 

Similar expression is obtained in the case of second filter 

y(i,j)=0,5z(i,j)-z(i-l,j) + 0,5y(i-l,j) (3.47) •. 

In using equation (3,46) and eqn, (3,47) it is to be 
remembered that the input is periodic and system is causal, 

i • e , , 


and 


x(i,j) = 0, Y'i < 1 0 


y(i, j) = 0, Yi < 0, j < 0 


For computing the impulse response array, we excited 
the system by a 2D delta function ^(i,j)* For the compu- 
tation of 2D Z-transform (on the unit circle) of second 


4 
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and third order noments of the output, we have used fast 
fourier transform routine FFTP (IMSL)» Fig* 3,10 through 
3,13 depict the 2D cyclic autocorrelation of the input, the 
actual and estimated amplitude and phase . response respe- 
ctively of the all pass system used in the study. In such 
cases we cannot take input noise array of arbitrary dimension. 
In order to maintain the idealness of autocorrelation, we 
have to arrange the elements of the sequence in a particular 
way to get the array. Considering all these facts, it seems 
that for such simulation purposes, computer generated random 
number are more convenient than shift register sequences. 

Such pseudo random sequences/arrays, can however still 
bo used as a substitute for computer generated random 
number for many applications. The periodic structure and 
almost ideal autocorrelation of pseudo random numbers help 
us in this respect. 

Since input noise is a pseudo-random array, the 
generated texture should be named 'pseudo random* texture. 






3,11: The Actual Aaplitudo response of the filter 
used in Sirnuiation, 





CHAPTER 4 


ANALYSIS AND SYNTHESIS OF TEXTURES BY fClNC/JUaX RF iODELS 
4.1 INTRODUCTION 

In the last two chapters, we have mentioned sample 
domain and frequency domain characterization of textures. 

In this chapter, vi?e will study a nev/ class of RF irodels 
where the concept of neighbourhood is a general one, i,e,, 
it can also take of non-causal neighbourhood. These models 
may be called as non-causal spatial interaction model as 
described by Kashyap [lo ]. Such RF models represent an 
image. as a distribution of gray levels over a finite or 
infinite lattice where each point is supposed to be equi- 
distant. 

Every such model is characterized by a set of 
neighbours, the corresponding coefficients and an independent 
identically distributed random sequence with prespecified 
probability density. Typically an image is characterized 
by 2D scalar data, the gray level variations defined over 
a rectangular or square lattice. In such lattice structures, 
gray level at a particular lattice point has a statistical 
dependency on the points inside a specified neighbourhood. 

There are two nonequivalcnt ways of specifying the 
interaction between pixels in a given lattice, leading to 
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two different classes of models, the simultaneous models 
[22 3 and conditional Markov models. Here we will 
discuss only the spatial interaction models over square 
lattice. These models do not require that the underlying 
field is Gaussian, These models are therefore quite general, 
i*e*, a variety of image types can be generated by choosing 
different neighbour sets, IID sequences of different 
probability density functions and different parameter sets, 

4,2 INDTATIONS 

r = ^ ® “ (^ 2 ^, 32 ): Co-ordinates of grid points, 

r^^ and s^^ integers, y(s).= gray level or intensity at pixel 
position s, y(s) is a non-negetive, real number, 

_0_ ~ ~ ^ “ ^1*^2 ~ ^ 

= set of grid points of a finite image of size MxM. 
o 

I =s s = (s ,Srt), — <» < S, ,S« < “ 
y = col [y(l,l), y{l,2),,.,, y(M,M)] 

= vector of a finite MxM image 
N = a neighbour set, a subset of (0,0) 
z = ® pair of operators, defined below 

If r,s £ 1 ^, z^y{s) ^y( 

= a pair of frequencies in both directions. 
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^ i = frequency in ith direction 
exp(-f-ip\) _ (expV'-i;^^, X 2 ^ 


A(2) = 1 . 


x4n 


© z-* 


AQ(r)= A(z = e^27ir/M ^ 

Y(.) = d.f.t. of l^y(s), s ^ J7_^ 

U(.) = D.F.T, of s 

R(.) = Autocovariance function 
S(A) = Spectral density 
P(.) = Probability density 

The convention for labelling members of N are as follows: 

x(-l,0} North 

x(0,-l) West . x(0,l) East 

x(l,0) South 

For example if the set N is represented by the 
coordinate as given 

N= (-1,-1), (-1,0), (1,0), (1,-1), (1,0), (1,1) 


f 

afc 2 

1 

k— 3 

r 

4 



r* 

i 

r 

\ 

f. . ^ n 

'r.. 

1 

^ * 

V - i 

r — 


then pictorially. 
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4.3 SUMMARY OF DISCRETE SPATIAL INTERACTION MODELS 


In this chapter we will discuss only discrete spatial 
interaction models with gray levels either real values or 
integers. In all such models the input is a set of indepen- 
dent random variables u(s) with a common density|%^. Such 
sequences are called independent identically distributed 
sequence (IID). The random field (RF) image y(.) is 
generated by transforming the input sequence in a suitable 
manner. In general, the system representation of such 
transformation may be of the following type « 


u(.) 





y( .) 


Fig. 4.1: A generalized model for RF Textures 

Firstly, we will discuss about T 2 and then with proper 
context we will discuss Tj^. The system T 2 is characterized 
by transfer function A(z) which is defined by a neighbour 
set N and postulated to be 

A(z) = 1 - 2 © , 9 real (4.1) 

r^N ^ 

The system Tj^ and T 2 may bo Interpreted as a difference 
equation in sample domain or multiplication in transform 
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domain. The output of the first system corre- 

lated or uncorrelated depending on the system T^. Now if 
we want to fit such a model for image data, then our main 
objective is to estimate the finite dimensional parameters 
the probability of the input IID sequence. 

We can divide the generalized model (Fig, 4,1) in two ways: 

Case i t The infinite lattice random field (ILRF) models, 
where y{s) is defined over an infinitely extended grid or 
lattice, i,e,, for such image Y(i,j) - » < (i,j) < <» , 

Case 2 ; The finite lattice random field (FLRF) model is that 
one where y(s) is defined only over a finite region and in 
that case_Q 

In general for RF models whether it is a infinite 
lattice or a finite lattice we can represent the image 
data y as 


1 < i»j 


representing 


image. 


y(s) - a 


\ [y(s+ - a] + Yyu(s) (4.2) 

Ic 

s = (i,j) 

i j j rr 0 I Hh 1 f Hh 2 y # 


Depending on the statistical nature of u(s) we may obtain 
various types of RF models such as SAR (Simultaneous Auto- 
Regressive), OA (Conditional Markov) etc. One advantage of 
such RF modes is that they incorporate geometrical and 
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statistical dependency of the pixel on its specified 
neighbourhood, not considered in many other models. By 
such models we can generate synthetic textures obeying 
a given model with known parameters by using a set of 
random numbers. 

4.4 INFINITE LATTICE ivDDEL 

Let us consider Fig. 4.1 as the general model. In 
IL models the transformation T 2 is interpreted as a difference 
equation whose coefficients are from A (2) and 

y(s) - Z y(s+r} = l/p\5(s) (4.3) 

This equation may also be written as 

A(z) y(s) = V^a3(s), z = (2j^,Z2) (4.4) 

where are regarded as delay operators in the 2 directions 
and A(z) is given in equation (4,1). The different choices 
for Tj^ lead to different types of models. 

Simultaneous AR (SAR) Models ; 

For a SAR model defined over infinite lattice, 

we have 

T^ = 1 and v(s) = u(s) 

Simultaneous /IIMA (SARMA) Models ; 

For a SARMA model defined over infinite lattice. 


we have, 
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= B(z) 

= 1 + S f 0^^ are real 

rfN 


r - ' '^"r 


Equ ivalently. 


y(s) - ®r = ■^'\5(s) 


= VTp[u(s) + £ d u(s+r) 3 

T r N ^ 


Separable i'Apdels: 


Here A(z) and B(z) possess a factorization 


A(z) = A^(z^)A 2 (z 2 ^ 

B(z) = B^(z^)B 2 (z 2 ) 

The neighbour sets characterizing A(z) and B(z) are located 
in a quadrant. 

Conditional Markov Model ; 

In conditional Markov {QA) models, the neighbour set 
N, used in characterizing the polynomial A and B is symmetric, 
i,e,, rC;N=^-r^N, Secondly the system is implicitly 
defined by the following relation. 

S^(}t) = Spectral density of v(.) 

= lS[expV-l 2it)^, exp f-l 211^211'^ 


( 4 . 5 ) 



Equivalently, the sequence v(,) has the following 
covariance function, 

Ry(s) = 1 if s = (0,0) 

= -9^ if s^N ('^.6) 

= 0 otherwise 

For such iTOdels, the coeff icients Og^ follows the 
relation 

s 

VVe can assume the relationship between u(.) and v(,) as 
follows 

= VI(I) (4.7) 

To validate the assumption we can consider two mutually 
exclusive cases. 

Case i : There exists a polynomial A^Cz) so that 
A(z) = Aj^Cz^) A^(z“^) 

and 

A,(z) = 1 - 2 ©i z^ (4,8) 

In such a case Tj_ lin eqn, (4.7) is interpreted as a 

difference operator 
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So in that case, the conditional model reduces to a simult- 
aneous model. 


£ase__ia: Let us assume that ACz) has no factorization, so 
in that case eqn. (4,7) has no meaning. 

It is also important to note that a synsnetric poly- 
nomial ACz^) in a single variable z^^ is alw/ays factorable. 
The general nonfactoribility of a symmetric 2D polynomial 
A(zi,Z 2 ) is one of the key differences between the ID and 
2D theory, 

4.5 FINITE LATTICE MODELS 


These models were introduced by Kashyap [32 ] and 
they have many elegant properties. In FL models both the 
inputs and output (image) are assumed to be defined over 
a h\yM grid. 

As we are interested in finite lattice image so in 
the subsequent sections we will discuss only FL model, 
however, an equivalence of these with IL models will also 
be mentioned. 


Neighbour Sets; Markov Definition and Static narity 
Neighbour Sets : A neighbour set N is any finite subset 
of the total set ^1*^2^ — ^1*^2 — * Fer all such 

sets (0,0) ^N, A neighbour set is called symmetric if 


s 
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The first order and second order neighbour sets 
are defined below, 

^ XXX 

X • X . X 

^ XXX 

N: First order N: Second order 

In the last two chapters we have analysed texture 
where neighbourhood was assumed to be causal* A causal 
neighbourhood has all its members in a quadrant 

N = (0,0) N, i < 0, j < 0 

X x 

X , X X . 

X , X 

N: Causal N: Serai causal 

N is called as one-sided or unilateral if N is a finite 
subset of a half plane S"^ defined as 

Definition : (Half Plane S"^): is recursively defined 

as follows 

i) f ^2^^^ ^ ^1 
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From the above definitions it is seen that is not 
unique. Some examples of are 



The non-causal neighbourhood is a more general neighbourhood 
and their use is more realistic for discrete image data 
modelling. 

Markov Definitions ; 

ID Markov Sequences ; 

A unidimonsional sequence ^y(t) , t= 1,2,...^ is 
said to be markov of order m if 

p(Y(t) iy(t^),^<t)=p(Y(t)|y(t-l),...y(t-m)) ( 4 . 10 ) 

probability 

vjhere p is the £ density function. Here m is 

called the memory of the process. A ID process y(t) is 
called bilateral markov of order m if 

p(y(t)laii y(tj.), tj, t) 

= p(y( t) }y(t+l), y(t + 2), ..«,y(t+m)) 

( 4 . 11 ) 
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In a similar manner, we can define vector markov 
sequence. In such a case, every component of the vector 
is said to be projection of a markov process, 

2D Markov Definitions 1*351 

A 2D sequence y(s) is said to be bilateral markov 
with respect to a symmetric neighbourhood N if 

p(y(s)|all y(s*), s 4 ^ s‘ ) = p(y(s) jy(s+s* ) j all 

(4,12) 

Comparing eqn. (4,11) and eqn. (4,12) we see that 
bilateral markov in 2D is just a generalization of bilateral 
markov in ID and in general bilateral markov in 2D is 
known as 2D markov sequence. 

Stationarity and Isotropy ; 

A sequence y(s) is covariance stationary if 

Covariance [y(s), y(s+r)3 = function of r only 
^R(r) = R(-r) 

where covariance [xy]A E[(X-E(x))(Y-E(y))] 

A sequence y(,) is weak stationary if it is 
covariance stationary and 

E[y(s)]=^ys 

The sequence y(.) is called a strong stationary sequence if 
its probability density function obeys the following 


conditions. 
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p(y(s«) = x^, y{s»») = x^) = fCx^jX^^s^-s* • ) 

In this study by stationarity we will mean the 

weak stationarity in the covariance sense, A stationary 
sequence will be called isotropic if, 

R(s = (i,j)) = RCliUljj) = R(j,i) 

For an infinite lattice covariance stationary R.F, the 
spectral density S(^ is defined as 

S(}y) = Z R(s) expCV-l 211 s, 

2 — — 
sti 

t>d cC 

= s s R(Sj_,s )exp[f-l2?^(s 

4.6 ANALYSIS OF FINITE LATTICE MODELS 

Vie have discussed in the previous sections some of 
the properties of infinite lattice RF(ILRF} models. The 
greatest impediment to the use of RF models defined over 
infinite lattices for analysing finite images is the 
extensive computation involved. Even the computation of 

9 o 

some of the basic quantities like M xM covariance matrix 
for an MxM image requires extensive conputations. Conse- 
quently, most of the RF models, used in image processing 
use simplistic assumptions like isotropy, causal neighbour- 
hood etc., so as to reduce the computational complexity. 
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Of course, there is considerable interest in the 
approximation methods. However, every approximation method 
seems to be tailored to a particular RF cod del and a 
particular application. 

The approximation of a covariance matrix by a 
circulant matrix or by block-circulant matrix is justified 
only by asymptotic considerations, i.e,, the size of the 
image tends to infinity. These approximations are not 
suitable for synthesis of any image, obeying a given RF 
model. 

Here in this section, we will analyse the basic 
features of finite or lattice RF (FLRF) models for images. 

We will discuss the development of several familieis of 
RF models defined only over a finite square lattice, say 
which is Same as that of the given MxM image to be 
analysed. The basic idea behind the FLRF nK>del is as 
follows: 

Using a given neighbourhood set N, characterizing the 
model or system, partition the given image (MxKi) lattice 
into two mutually exclusive subsets XI j and Xlg* The 
boundary set XI g has every grid point which has at least 
one neighbour defined by N not in J2- . is the interior 

set, equals In FLRF models the eqn. for generation 
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of y(s), and y(s), a^e different. %)ecifi- 

cally every pixel in_Q_j follows the eqn. (4,2) i.e. 
they are same as that of ILRF model but each pixel in _n. g 
follows a modified equation as shown by eqn. (4,13) 

y(s) = - £ + V?u(s) (4.13) 

Thus the equation for the image vector y can be written as 


B(e) y = ii (4.14) 

where B(©) is M^xM^ and u is a M^xl vector of primitive 
random numbers independent or dependent on the basis of 
SAR and CM models. By choosing app^piately, 

the matrix B(0) can have eigenvectors, either fourier 
vectors or discrete sine or cosine vectors as the case 
may be. The distinction betweenJV^ andjQ_g can be given 
as follows 


SI- 


B 


s=(i>j)» s^A-and is+Sy^)-£ St for at least 
one number Sj^£:N J 


Sir - js=(i,j),s^Jl and (s+s^) ^ 


The eqn. (4.13) was chosen for the pixels y(s), 


such that 
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x*et they do not involve 

any where s» 

the overall set of m 2 equations are easy to analyse. 
SAR Model 

On the basis of the above mentioned partition of the 

fxnrte lattice Jlmto andJI ^ we can write the eqn. 
(4.2) as 


for s£ i\^ y(s) obeys (4.2) 
but for s£JVg 

y(s) - a = - e^[y^ (s + k) - q]+V3u{s) ( 4 . 15 J 

whe re 

y^Cs+k) with s=(i,j) and k=(p,q) 

= y{s+(p,q)) if (s+(|>,q)) ^ JV 

= y[(p+i-l)j^+l,(q+j-l)j^+l] (4.16) 

if s+(p,q) ^ Jl_ 

This modified eqn, (4,16) is based on the assumption that 
the lattice is a torus one, i.e,, for such lattices the 
opposite edges are similar. This is a very useful conception 
used in this and subsequent chapters, • 

Eqns. (4,16) and (4,2) give eqns, relating the 
image variable y(.) and the IID variable u(.). By using 
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the Same relation as given in eqn* (4*14) we see that for 
toroidal assurtption the matrix B(9) becomes a block cyclic 
of the form. 


B(©) := 


®11 ®12 ••• ®1M 

Bii ... 


B 


12 


# • « « • < 


B 


11 


where each B^^. is a cyclic submatrix and depends on the 
parameters For example, when the neighbour set of 

dependence N is 

B 

B 

B 


-1,0), 

(©,1), (1,0), (0, 

-i)j 

11 ^ 

eye 


m • m f ^ 

12 

eye 

^'^1,0 * 

0 3 

‘iM ~ 

eye 


.. 0] 

i. . = 0 
ij 

1 for 

j 1,2,M. 



To ensure stationarity, the coefficients must obey 

^s(a) 4 (1- (4.17) 

where, 

= col [exp[f-l ^ {s-l)’^r, r ^N]? 1 =(l,l) 
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In a similar manner, we can construct the simultaneous ABMA 
or SARMA model. 

Conditional Markov i(M) Model ; 

In previous chapter we discussed the properties of 
a CM model defined over an infinite lattice. Again, we can 
have CM models defined over a finite lattice. Vie can 
also tackle the problem of edges by toroidal assun^tions. 
The vector matrix relation in CM model is 

H(©)i = e (4.18) 

where H(©) is a block cyclic and symmetric, ^ere to 
ascertain stationarity we need 

(1-2©.'^ efg) >0 V* sellu 

© = col [©^, r^N] 

and 

= col [cos |2S(s-l)^r, r ^ Ng] 
where N is the asymmetrical half of N,i*e,, 

N =: Ng U 

Ng=|^3: : -r^Ng j 

NgnUg-o ■ , 



85 


4.7 properties OF FL AODELS 

Considering the matrix B(9) defined in eqn. (4,14) 
we can derive some informatiore about FL SAE model from its 
toroidal property. 


l) Since B(0) is a block cyclic matrix ^ hence the eigen- 
vectors p.^.(0) of B(©) are the fourier vectors f - where 

J ij 

***/^i ^ ^ (Mxl vector) 

and 

~ Xj ] ( Mxl vector) 

with == exp [ipi ^ ] 

So given the set of parameters 
the eqn, for p. .(©) as 

XJ 

ny(e) ) 



we can write 




(4.19) 


where k = (p»q) 


2) Again from the vector-matrix relation we can write the 
covariance Q of Y is as 


Q= cov [y] = [B^(0) B(e) 
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R(k|i^) — autocorrelation at lag Ck,4.) is given by 


R(k,e) E )/ll>i.(©)ll^ (4.20) 

lfj=l ^ J AJ 

The discrete spectral density as defined by the DFT (2D) of 
R(k»^) and then for such block cyclic matrix B(6) we have 

=f/|| J e^^2)||2 

= _L 


for iP =: 


3-J 


'^=222 
^2 M 


CM Models; 


Since in FL toroidal Q*A model H(9) is a block-cyclic 
symmetric matrix, some properties of CM(FL) models are 

1) Similar to FLSAR model, the matrix H(0) has fourier 
vectors as its eigenvectors. If we assume that (©) is 
the eigenvalue of H(0) corresponding to f^j we have 

2) Without any loss of generality we can assume a =: 0. 

Then, y= [HCo)]"^ (4.22) 
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where e is a correlated sequence. 

To obtain the covariance of vjq multipily eqn, (4,22) by 
and take expectations on both sides, 

Q = [H(9)3“^ E[y^e] ^ 

-I 

= [H(9)] y* (4,23) 

from eqn, (4,6). 

and cov[e3 = E [ee"^] = H(0) (4,24). 

4,8 RELATIONSHIP BETWEEN SAR AND CM .MODELS 

f 

In general, two FLRF models defined over same 
lattice _|7- said to be equivalent if their corresponding 
output vectors are identical. Unless otherwise stated# 
equivalence means equivalence in all respects. 

Direct comparison of SAR and CM models are hampered 
by the fact that input e to the CM system is correlated. It 
would be desirable to modify CM model to a new model where 
the input is uncorrelated like that of SAR* This new 
representation is also helpful in synthesis of image 
following a given Cm model. 

The CM model as given by eqn. (4,22) has the following 
representation 


VHte) M 


(4.25) 
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IW) 


M 

I 




f 


ij 


w 


(4*26) 


where u is the . normal zero mean IID sequence and V*H(9) is 
real and positive definite. 


Eqn. (4*25) again may be written by using eqn, 
(4,26) as 


M 

1 , 


1*4 _ . 

i.ii 


(4.27) 


So e in eqn. (4.22) is related to u by the eqn. 


e = VHl©) u 


(4.28) 


Even though eqn. (4,25) has a form similar to that 
of a SAR model but, it is not necessarily a SAR model, A 
FL CM model represented as 

GCz^z^) [y(s)] = V^e(s) (4.29) 

and it is equivalent to a SAR system represented by eqn, 

B(©‘)x = Y^ii (4.30) 

if G( 2 j^,Z 2 ) passes a factorization given in eqn. (4,31) 

G(zijZ2) = 5 1 1 a(z3^,Z2)} 

= k A(z^,Z2)A(Z^’“^, z^**^ ) (4,31) 

If GCzj^fZ^) does not follow eqn. (4.31) but there exist 
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polynomials and and a constant k 

such that 

S(Zi.Z2>ll“'l(2i'Z2)!l^ = k||o2(zj,Z2)|l=^ (4.32) 

then the corresponding CM model has an equivalent SARMA 
model. If G(z^^Z 2 ) factorizable then the CM 

model has no equivalent simultaneous model, 

4,9 PARAMETER ESTIMATION 


Suppose we have already specified the model type and 
the neighbour set. Let 0 = (£> ^ ) be the vector of unknown 
parameters, where 0 is the vector of parameters and 
is the noise variance. The problem of estimating ^ from 
the given image vector y is primarily a statistical 
problem. This problem can be solved either by least 
square (LS) method or maximum likelihood (ML) method* In 
the case of LS estimate the estimates are given by 


0 = [S z(s) z^Cs)]"^ (r z(s)y(s) 

and 



(4.33) 


(4.34) 


where 

z(s) = col Cy(s+r), r6N] (4.35) 

It has been shown by Ord [37] that this method does not 
give consistent result for nonunilateral neighbourhood. 
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MaximiM Likelihood Estimate (ML£) i 


In case of Fl images, the WLE gives a consistent 
and efficient estimates of the parameters. For such MLE 
we have to find the log likelihood function of the output 
vector Y.» To obtain such an expression we impose the 
condition that input IID sequence is a zero-mean Gaussian 
which gives a Gaussian output as the system is linear. 

If we represent the system by its tranfer function A( 22 ^^ 22 ) 
or by its matrix from B(©) then likelihood function comes 
out to be 


P(l I © J exp[- 2 A(z 2 ^,Z 2 )y(s)]^] 


( 4 . 36 ) 


where J is the Jacobian 
or. 


In pCxI^y'f ) - B(©)] - ^ In (2ii ^ 

- S (y(s)-9^z(s))^ 

s ^ A. 

since det [B{9)3 = (T“ ^ 

we have 

2 

In pix I ^ ^ Tn(2iiP) 

^ 2 {y(s)- 0^ z(s))^ 


( 4 . 37 ) 


( 4 . 38 ) 
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The ML estimates are obtained by maximising eqn. 
(4,35) w,r,t, 0 and -p* • Since log likelihood function is 
a quadratic in 9, the estimation may be done by numerical 
optimization schemes like Newton-Raphson etc., which are 
computationally expensive * 

Here we shall describe an iterative schen^ for 

estimation which is an approximate version of the MLE. In 

such a scheme we approximate expressions like ln(l-i*a) by 
2 , 

a ~ a /2, The corresponding log likelihood function 
denoted by J(0, P) becomes 



J>) = -v\ + 0.5 - j- In {2»l'p) 


where, 

- 5 ^ E Cy(s) - 

(4.39) 


V = S C , m X 1 vector 

(4.40) 


R = 2 (S sj - S ^ matrix 

(4.41) 


Cg= col [cos ^ 

(4,42) 


Sg= col [sin ^ 

(4.43) 

After 

defining the matrices and vectors the 

algoritlim 


for approximate MLE is as follows: 
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®t+i = (V - ) 

and 

ft = ^ ^ (y(s)- ejz(s))2 

where 

S = E z(S) z(s)^ 

(mxm matrix) 


U = E z(s)y(s) 
s {IJI^ 

(mxl vector) 


( 4 . 44 ) 


( 4 , 45 ) 


( 4 . 46 ) 


( 4 . 47 ) 


The initial value 0^ is chosen as U, and here m 

is the dimension of 0, 


CM Model ; 

An estimate with good consistency and efficiency may 
be obtained in CM models by ML estimates. But it is not 
easy to derive an explicit expression for the log likeli~ 
hood function for such CM models because Jacobian of the 
transformation matrix is difficult to be evaluated. By 
assuming a toroidal structure and Gaussian correlated 
sequence e(s), the log-likelihood function may be given as 
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in p{x|a.f)= j|ln (l-2el^ ).^ln(2.f) 


- ly 5C.^H(e)y 

(4.48) 

We see that the exponential term in In pCyj©,')^) is linear 
9l unlike that of SAR model, A computationally efficient 
and consistent estimation procedure is as follows 


* r 

© = [ 


S q(s)q(s)^]“^ ( £ q(s)y(s)) (4.49) 

s J s ^ Jft 


and 


y* = ^ S (y(s) ~ 0 * T q(s))2 

M s4ji.^ 


(4.50) 


where JD-j is the interior region. This procedure is 
suggested by Woods [36]. . In our study, estimation 

of ■ parameters for both SAR and OA followed the above 
mentioned scheme and experimental results show that the 
estimates are quite close to the actual values and moreover 
this scheme works well for non Gaussian cases also (see 
section 4,11) . 

4.10 CHOICE OF NEIGHBOURHOOD FOR SAR AND CM iVDDELS 

From ID time series modelling, it is known that 

order 

only a model of appropiate £ gives the best result for a 
set of discrete data for forecasting and prediction. A 
similar situation is true in 2D spatial processes, but, the 
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problem with 2D data set is the variety of models available. 
So if we want to reconstruct any image from the parameters 
extracted from it, then the quality of reconstructed image 
depends on how the underlying model is similar to the true 
model and here comes the importance of the choice of 
appropiate neighbourhood. 

Though, one way of picking up the appropiate model 
is by visual inspection, we give below a quantitative 
measure for this. One way is to do pairwise hypothesis 
testing. But such acheme is not consistent and also not 
transitive. If a model is preferred to C2 and is 
preferred to then it does not follow is preferred to 
• Such a decision rule is not consistent in the 

sense that the probability of choosing an incorrect model 
does not go to zero if we take infinite number of 
observations. The other choice is based on Bayesian 
decision rule [35]. 

Bayes Rule for Choosing SAR Model : 

Suppose that we have three sets and N3 of 

neighbours having m^jm^ and members respectively. 

For each model we can assign one SAR model 

y(s) = Z ©.y(s+r) + u(s) s^Jl^ 

y(s) = S \ yj(s+r)+irpj^u(s) s^Jl^ 
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where the significance of such equations are already 
defined. 

The models C^, i = 1,2,3, are mutually exclusive. 

By Bayesian rule we can say that the optimal decision rule 
for minimising the average probability of decision error 
is to choose the model that maximises the posterior 
probability P(Cj^ly) , where ^ is the observation vector. 
The quantity P(Cj^jy) is computed by Bayesian rule, 

j^)P(C^)/p(y_) , We can set P(Cj|^) same for 
all i in the absence of any contrary information, so 'Uiat 

Pil |Ci,) =/p(y|©,j^) p(0, j>lCj^)id©i d|> 

This expression, with the Gaussian assumption is similar to 
the expression which is obtained in the section on para~ 
meter estimation. By assuming that p(9,'p/C^) is a-priori 
density then by using asymptotic information it can be 
shown that 


P( I 




M ^/2 


In P 


. I 
^ 2 


S 


ln(l-©“% 


ks-^^k 


-T 


^ks 


®k) 


m. ^ A . X 

- — “ ln(M )+ In Y ^ * 

In general, computation of such an expression for ail values 
of k is not computationally efficient so here we will give 
an approximate test statistic by dropping the terms due to 
prior densities. Then the decision rule for the choice of 
appropriate neighbour set is following. 
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where 


Choose the neighbour set N if 

k 

k = 1 min 5 C ^ j 


'k 'j ~ ^ ln(l-2 + e“^ §^)+ 

+ nij^ In(M^) 


( 4 . 51 ) 


Qks= 


S |^3 = col [sin 2a ((s-l)Tr), r^fj] 
Bayes Rule for Choosing CM Model ; 


Suppose, we have three sets Nj^»N 2 ^3 neigh- 
bours having 2m^, 2m^ and 
ctively. For any nwdel Nj^, we can write. 


2 m 2 members (symmetric) respe- 


y(s) = E y(s+r) + l/pj^ e(s), s £Sl.^ 

^ 6 ,^k 

= Skr f k =^-*^6 

r i 

where the symbols are defined previously, and e(s) is 
correlated Gaussian sequence. Then in the similar way^ 
of SAR> we can say that the approximate test statistics 


will be: 
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Choose the set N. if 

iC 


“ = I T i c. 


where. 


Ck = 




+ In (M^) 


where , 

% = col [e. . k eN ] 


and 


1 r ^ ^ r' T N 7 

'^ks 


col [cos ^ (s-l)^r), r^N, 3 


(4.52) 


where is the asymmetrical half of I^, i.e 


• » 


N, = N U R 

I' \ =k 


N, 


r : -r N , N 0 ^ 

®k \ ®k ®k 


In the next section where we will discuss the simulation 
results, we will show how test statistics detects the 
actual model for synthetic images for both SAR and CM models. 


4.11 EXPERIMENTAL RESULTS 

To see how the SAR and CM models and associated 
estimation scheme work, different images, using different 
models have been generated and displayed for visual tests. 
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The fast schemes for generating toroidal synthetic images 
by a known SAR and cal models are given below» 

Generation scheme for MxM imAge (SAR) ; 

i) Generate a sequence of IID numbers (Gaussian or 
uniform )of length M and arrange then to form a MxifA 
array. Let that array be U . 

ii) Given a set of neighbours find the 

eigenvalue matrix p, of the system as 


Pi3.(0) 


1 + 


Z 

r eN 


0 z. 

r 1 


t 


where r= 




iii) Take the 2D D.F.T. of the array U. Let the 
transformed array be 

iv) Generate an array Ey where 

Fy(i»j)= (^u (4,53) 

v) The final image obeying the given SAR model is 


y = 



F"^(Fy) 


= ^ j |a( j ) ) JP 3 

M ^ 

In actual practice we did all the transforms by the BISL 
routine FFT2 and FFRDR2* 

The above algorithm for generating image was used 
with different neighbourhood (pre-specified) and with IID 
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sequences with different densities. Parameters were 
estimated from the generated image by using approximate 
ML estimates, in any random number generation, we require 
some initial values (called ISEED in IMSL) and it has been 
found by experiments that estimated parameters change with 
different ISEED value and that change is also quite 
random. 

Generation Scheme for h\m Image fO/iV ; 

Since the generating eqn, in CM model can be 
written as 

rH(9) X 

hence the generating strategy for CM is same as that of SAR 
except the eqn, (4,53) where in CM, eqn. (4,53) is nidified 
to 

Fy =(Fu (i,j)/ Vp(i,j))sl^ (4.54) 

Hence also we take uniform as well as Gaussian density 
for the IID sequence , 

Results for SAR and OA model parameter estimations 
for different ISEED value will be given in tabular form. 

In the next set of tables we will give the estimated 
parameters with the test statistics for finding the 
optimal neighbourhood. Here we have considered only 
noncausal neighbourhood with 2,4#^ and 8 parameters for 
both SAR and CM models. 
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GAUSSIAN SAR ESTIMATED 
ACTUAL (6 Point) 


r 

= 1*1111 

0.9982 

1.035530 

1.081630 





®-ii 

0.28 

0.2988 

0.3012 

0.2877 

© T T 

Q 


d 






- 1-1 

-10 


^10 

0.12 

0.1372 

0.1464 

0.1220 

0 

> 

► . 

a 

®-l-l 

- 0.14 

0.1612 - 0.1589 

- 0,1729 

a 

> 

< 


o 

®11 

- 0.14 

- 0.1612 - 0.1589 

- 0.1729 



.0 

®11 

®10 

0.12 

0.1372 

0.1464 

0.1220 





® 1-1 

0.28 

0.2988 

0.3012 

0.2877 






GAUSSIAInI SAR (4 POINT ) 


^10 

®0t1 

®01 

®10 


= 1.1111 

1.02123 

1.07213 

1.00213 

0,36 

0.4123 

0.4005 

0.3823 

- 0.12 

- 0.1423 

- 0,1252 

- 0.1512 

- 0.12 

- 0.1423 

- 0.1252 

- 0.1512 

0.36 

- 0.1423 

0.4005 

0.3823 


0 


-10 




0-1 


®01 


& 


10 
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UNIFORM (SAR) 6 POINT 


j ) = l.llll 

1.021346 

0.998312 

1.021672 

w 


0.28 

0.3087 

0.3054 

0.3124 



0.12 

0.1423 

0.1482 

0.1480 , , 

© 

© 




- 1-1 

-10 

-11 

- 0.14 

- 0.1642 

- 0.1601 

- 0.1650 



• 

O 

I 

- 0.1642 

- 0.1601 

- 0.1650 



0.12 

0.1423 

0.1482 

0.1430 © 1.1 

~10 

®11 

0.28 

0.3087 

0.3854 

0.3124 





CM (6 POINT ) GAUSS 



= 1.1111 

1.2312 

L. 01245 

0.99824 


.2800 

0.2912 

0.2992 

0.3012 

»-10 

®-10 

0.1200 

0.1421 

0.1291 

0.1282 

» 0-l ®01 

^— 1.1 

— 0 . 1400 

- 0.1232 

- 0.1501 

- 0.1542 


®11 

- 0.1400 

- 0.1232 

- 0.1501 

- 0.1542 


®10 

0.1200 

0.1421 

0.1291 

0.1282 


©V 

0.2800 

0.2912 

0.2992 

0.3012 
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CM (4 POINT) GAUSS 


ACTUAL ESTIMATED 


®-l-l 

0.28 

0.2942 

0.2852 

0.2915 



®-ll 

**0 # 14 

-0.1240 

-0.1493 

-0,1413 


®-ll 


“0. 14 

-0.1240 

-0.1493 

-0.1413 




0.28 

0.2942 

0.2852 

0.2915 



f 

1.1111 

1.0243 

1.0823 

1.1056 

®1-1 




CM (4 POINT ) UNIFORM 



0.28 

0.2752 

0.2901 

0.2713 

®-l-l 

®-ll 

-0.14 

-0.1482 

-0.1391 

-0.1542 



— 0 . 14 

-0.1482 

-0.13'91 

-0.1542 



0.28 

0.2752 

0.2901 

0.2713 


®11 

1.1111 

1.1623 

1.1192 

1.2109 
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Computation of Parameter for an Image (following SAR 
Model) when the real neighbourhood is N = (0»1), (0,-1), 

(1,0), (-l,0).0riginal ©j_q=0_^^^=O.365 -0.12j 

I = l.llll 


& 




Positions 


1.0355 


0.4858 

0.4858 


179.42 


(- 1 , 0 ) 




( 1 , 0 ) 


1 .0052 


0.4005 

-0.1152 

-0.1152 

0.4005 


. (- 1 , 0 ) 

80.76* 

. ( 1 , 0 ) 


( 0 , 1 ) 


1.0686 


0.5616x10 
0 .4041 
0.5550x10 
0.6550x10 
0 .4041 
0 .5616x10" 


-1 


(- 1 ,- 1 ) (.. 1 , 0 ) (- 1 , 1 ) 


-1 


1 

--- j 

► — t 



5 G 

* 


/ 

Ik... . „ . 

< j 



|r 

. . jj 



587.48 


( 1 ,- 1 ) ( 1 , 0 ) ( 1 , 1 ) 


1.0143 


0.6963x10 
0 .3966 
0.7091x10 
-0.1067 
-0.1067 
0 .’7091x10 
0.3966 


-2 


-2 


230.07 


( 0 ,- 1 ) 

( 1 ,- 1 ) 






^ ^ 

"‘“Q 

j 


w 


i' - J 


lb«»> 


} ■ £] 

/ - 

'"*1 



( 0 , 1 ) 


( 1 , 0 ) ( 1 , 1 ) 


-2 


-2 
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Computation of the parameter Cj^ for an image (following CM 
model) when the real neighbourhood is N = (1,1), (-1,-1), 


(-1,1). 

Original ® ]_j_ = 


p = 1 

.1111 


®ii = 

8 - 1-1 = 0-28 




1 

0)1 

o 

.‘k 

Positions 

1.1292 

^=0.3221 

420.72 

- 1-1 

11 

1.4213 

e _^ ll =042121 

621.63 

1-1 

-11 


© .1512 



-11 

1.0321 

0^^ =0.2922 

121.43 

1-1 

11 


0 j ^ q = O . O 212 



-10 

1.1021 

412.76 

0-1 

01 

@0 #0092 


10 
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Real Image Synthesis using SAR Models ; 

Our discussions so far have been concerned with 
synthetic image patterns generated by known SAR models. 

We conclude this chapter by studying the appropiateness 
of the SAR models to real images. For this purpose, we have 
taken a 32x32 portion of Monalisa image. Since the entire 
image may not be homogeneous we have divided the entire 
image into 4 blocks of size 16x16, Parameters 
were estimated for each block by the approximate ML estimate 
scheme with 12 point neighbourhood. The neighbourhood 
includes the points (1,0), (1,1), (0,1), (-1,1), (-1,0), 

(-1,-1), (0,-1), (1,-1), (2,0), (0,2), (-2,0), (0,-2). We 
have also estimated the parameters considering 8 point 
neighbourhoods. In that case, the neighbourhood N contains 
( 1 , 0 ), ( 1 , 1 ), ( 0 , 1 ), (- 1 , 1 ), (- 1 , 0 ), (- 1 ,- 1 ), ( 0 ,- 1 ), ( 1 ,- 1 ). 

In the case of 8 point neighbourhood the quality of the synthesize 
synthesized image was not good. From the estimated para- 
meters and sample means, we compute the residuals m as 


w 


^ [H(0) 1 - a 1] 

T 

where y. ’the block image and a is the sample mean* From 
this sequence w we have generated a set of random numbers 
whose histogram was approximately same as the histograms 
of the residuals. The new sequence and the parameter 
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set were then used to synthesize the image* The synthesized 
image closely resembles the actual one except the edge 
effect of the block operation* 

4.12 DISCUSSION 

In this chapter we have given a complete methodology 
for analysing a given set of spatial data. We have also 
discussed the estimation scheme for SAR and CM raadels from 
the computation of test statistics. We have not said 
any thing about the superiority of SAR over CM model or 
vice-versa as the appropiateness of one class of model 
is entirely dependent on the data itself. 






CHAPTER 5 


TEXTURE SEGMENTATION BY SAR MODEL 
5.1 PRINCIPLE OF SEGMENTATION 

Segmentation is a process that breaks up a conposite 
scene into its constituents. Though there are several 
approaches for tackling this problem [3 ] ys't it is an 
active field of research because of its importance as the 
first processing step in any practical scene analysis 
problem. In natural composite textures (e.g. combination 
of land and sea) we may or may not have any sharp boundary 
as the line of demarcation. In such cases v4iere the 
contour is missing, it has been observed that statistical 
measure gives an estimated contour which is very close to 
the original one. Here segmentation by gray level 
thresholding will not work satisfactorily because there 
may be two regions having the same brightness average 
but slightly different variances so that it is not distin- 
guishable. In a composite image (texture) having two 
homogeneous regions we assume that all the pixels in one 
region are having same statistical properties. In fact, a 
composite texture may have more than two homogeneous 
regions but solving a two region problem can always be 
generalized for n-region textures. 
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Though there are several statistical segmentation 
schemes, the proposed scheme for segmenting a random field 
composite texture based on spatial interaction model is 
superior to them in the sense that it is computationally 
very simple and to some extent accurate. In this chapter 
we will discuss the principle of such segmentation scheme 
and the associated results, 

5.2 APPLICATION OF RF NDDEL FOR SEGMENTATION 

Given a composite texture X, we assume that for each 
homogeneous part of it we can have a RF model and model 
parameters are different from one region to another. Let 


X = U (5.1) 

where and X 2 are the homogeneous random textures. Since 
and X^ are following a toroidal SAR model, we write 


X- I YiCs) + S 0, y-iCs+k,) - l/^u(s),s^/2.j 
^ 1 k,ON, H ^ ^ ^ . 

(5.2) 


' 1^1 


yi(s) 


+ E 


9, 



(s+k^)=Yp^u(s) st-h-g 


X^: 72(5)+ 0 j^ y 2 (s+k 2 )=rf 2 ''(s) 

^2^^2 ^ (5.3) 


similarly, 
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where the various parameters appearing in eqn. (5.2) and 
eqn. (5.3) are explained in the previous chapter. The 


entire image is composed of two regionsj region~l and 
region-2, represented by and 50^ ^ respecti- 




After estimating 


vely, where 9, =2 9, 

1 

the parameters for each region we use them to find the 


similarity of them inside a region. The segmentation 
procedure begins by breaking the entire image into smaller 
square blocks. We then consider each block image as an 
independent one and the parameters for each block is esti- 
mated, The Same procedure is followed for all such blocks. 
We thus obtain a set of parameter for each block. The 
final step is to cluster these set of data by any standard 
scheme. The clustering of these clata will determine the 
segmenting line. Though the entire image may be correlated, 
for the purpose of segmentation, we assume the block 
images to be independent. This is not a valid assumption 
and will give some error in the final result. But 
experience has shown that this error is negligible with 
moderate block size (8x8) and it varies inversely with the 
window size. The entire process can be visualized by 


Fig. 5.1. 

Here we consider a MxM image with windows of size 
NxN. So, effectively, we will have K such blocks where 



110 



segmenting line 

Left portion of line 
with (©^, -p^) and 
right with f 2^ 


Fig* 5,1: A Block Composite Image 


H 


K - M /N^. V/e estimate the parameters of each block, 
treating them as independent, and arrive at K sets 
Then by clustering these parameters we can separate the 
non-uniform blocks. Since clustering in 2 region problem 


is a binary decision making process, the estimated curve 
will always be a staircase one. Evidently, the estimated 
curve will be close to the actual one if the block size 
is reduced but, we will show that it has other implications 
also. This method is tested with different segmenting curves 
like straight lines, parabolas, circles and for each case 
satisfactory results are obtained. 


5,3 PRINCIPLE OF CLUSTERING 

In this section we will discuss the principle of 
clustering and the algorithm used. Clustering is also 
known as unsupexvlzed learning. Roughly speaking, clustering 
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procedures yield a data description in terms of clusters 
or groups of data points that posses strong internal 
similarities* There are several methods of clustering for 
a set of data but all are heuristic in the sense that they 
do not possess any mathematical properties but are easy 
to use [38 ] « Since the clustering procedures have no 
mathematical foundations, they depend on the nature of the 
data* For all algorithms at least some properties of the 
data should be known a priori* 

In this study, the algorithm used for the segmentation 
procedure, is known as K-mean algorithm [39]* In this 
algorithm we assume that the number of clusters at final 
stage is K. This is one approach of general iterative 
optimization scheme. We will describe the algorithm and 
also give an example of it. 

The K-mean algorithm is based on the minimization 
of performance index which is defined as the sum of square 
errors. This algorithm consists of the following steps. 

1, Choose K initial cluster centres z^(l), Z 2 (l),...»Zj^(l)* 
These are arbitrary and are usually selected as the first 

K samples of the given sample set. 

2. At the Kth iterative step, distribute the sanples-^^ ^ 
among the K cluster domain using the relation 
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X4Sj(K) if 

l|x - Zj(K) II < II X - Z^(K)1|, V'i=l, 2 ,...,K 

ij^K (5.4) 

whor© Sj(K) denotss th© set of samples whose cluster centre 
^j(^)»Ties in expression (5.4) are solved arbitrarily, 

3. Compute the new cluster centre Z.(K+l),-^ j=l,2,.. .,K 
from the result of step 2 such that the sum of square 
distance from all points in S.(k) to the new cluster centre 
is minimised. In other words, the new cluster centre 
Zj(K+ 1) is computed so as to minimised the performance index 
J j , where , 

J. A S 1 lX-Z^(K+l)! j=l,2,...K 

^ =^X€S.(K) ^ ^ 

J (5.5) 

The Z.(K+1) which minimises J. is simply the sample mean 
J J 

S-(K)* So, the new cluster centre is given by 

%J 

Z.(K+1) = ^ E , ^ » 3=1, 2, ...,K (5.6) 

J X Sj(K) 

where N. is the number of samples in S.(K). The name 
J . 

’K~mean* is obviously derived from the manner in which 

the centres are sequentially updated. 

4. If Z^(K+1) = Z.(K) for j = 1,2, ...,K, 

the algorithm has converged and the procedure has 
terminated. Otherwise go to step 2. 
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Convergence of such a scheme is not proved and 
convergence depends oh the value of K and the type of data, 

V^fe will also give an example of such scheme for 
a 2D data * 


Example : 

1. Let K=2 and let 
us take 

Zj^(l) = 

and 

Zryil) = [ 0 3=^2 ^ 






WMM. 


X 

— * 

If > 






Xi 

^ ■ 







i. i 
















? > 

!/o ' 










1 



-J 







IZ 

M 







,v 



1 








i 

1 

! 

1 





2. Since 

lIX^L-Z^CDlUIIXj^-ZyU'l'l ‘ ’Fig.V2:'2D” Clustering 

and 

1 |X3-Zj_(l)l 1<I ! 

So, 


32 ( 1 ) - 

3. Update the cluster centres 


1/v V > r o# oo 1 

X = 2(Xj^+X3) = L 0.50 J 


z.(2) ^ ^ 


1 x^St(1) 


^ t ^ “ 18 •*■^^20^" ^5.63 ^ 

^'^2 x6S2(1) 


ZA2) = ^ ^ ^ “ 18 
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4. Since Zj(2) 4 Zj(l) V-j, so go to (2) 

2, With new cluster centres we get, 

||X^ - Z^(2)||<i|Xj^- Z2(2)1|,1^^= 1,2,..., 8 

and 

l|X| -- Z2(2)i|<||X^- Z;^(2)1 |,v£= 9,..., 20 
Si(2) = [x, ... Xgj 
SjCz) = ^i^Xg ... yiX 

3, Update the cluster centre 
z^(3) = I (X^ t ... + Xg) = 

Zj(3) = ^Xg+ ... + 

4, Since Z..(3) Z.(2), j=l,2 so go to step 2. Step 2 

kJ ij 

yields same result as in last iteration so Zj_(4)=Z^(3) 
and Z (4) = Z„(3). Step 3 also yields the same results* 
Since Z,(4) = Z.(3), j=lf2 the algorithm has converged 

w 

and it hias yielded the cluster centres 

^1 ^2 1 ^ 7 , 33 ^ • 

This algorithm is very sensitive to the initial 
cluster centres but generally the first two data in the 
data set are taken as the cluster centres for the first 
iteration* 
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Thouqh the above example deals with clustering 
of 2 dimensional vector but in practice for image segmen- . 
tation we need clustering of feature vectors of dimension 
n+1 (for n~point neighbourhood). 

5.4 SIMJLATION RESULTS 

In this section we will discuss the simulation results 

obtained by tho above scheme. Various composite textures 

wor<‘ genera tt?d using different SAR models, m have taken 

only 64x64 image but various' window sizes were used for 

sGipenta tion purposes. An empirical relation between the 

window size and error is given where error is defined in 

the mean square sense. If ©jjj represents the parameter 

A 

vector for a fixed window w and 0^ be its estimate then 
the error is defined as 

1 6 I ■ 

/^ = 0 - 0 

or, 

c. ^ 

where N is the length of vector 0, 

In the following, we present the method adopted for the 
segmentation scheme and the method of error estimation for 
the synthetic composite texture. The generated image has 
two parts or regions as shown in the Fig. 5,3, 


N 


Irsl 1 1 


(5.7) 
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Fig* 5.3; Composite Texture 


This image v/as actually generated by superposition 
of two different images such that, in the composite image 
left portion of the line (at 45°) has parameter 
and the right portion is with parameter 
actual simulation = [0.12 0.12]^ and x “ 1*1111 

= [0.36 0.36]^ and ^2 ~ 1*9999 was used. We partitioned 
the entire 64x64 image into blocks of sizes 8x8 and 4x4, In 
case of 8x8 blocks we get 64 such blocks and 256 in case 
of 4x4 blocks. In Fig, 5.3 the line passes through the 


blocks 8,15,22,29,36,43,50 and 57 and they are really the 
composite blocks. Parameters vjexB estimated for each blocks 
and all the blocks on the left side of the line yields nearly 
same parameters as is to be expected and same is true for 
the blocks at the right side. After parameter estimation. 
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we have a set of 64 (for 8x8) vectors where each vector 
has three components when we estimate a 2 point problem. 
These vectors are called feature vectors and they are used 
for clustering. Since vie have considered symmetric neigh- 
bourhood, the 3 point vectors reduces to 2 point because of 
the symmetry. These vectors then can be plotted on a 2D 
plane. Graphically, the separations are distinguishable. 

As mentioned earlier, all clustering schemes are heuristic 
and they have inherent limitations in separating data points. 
In the present case also, we end up with a few misclassified 
samples but they are easily detectable, Vs/e show a segmen- 
tation result for the same texture for 2 point and 4 point 
parameter estimation in Figure 5,4, 




2 Point 4 Point 

Fig. 5.4: Segmenting line estimated by different SAR 
models. 



liS 


In Fig, 5,5, we have given the actual and estimated 
segmenting curves. In order to study the error in esti- 
mated parameter with window size,- v;e have taken a pure 
homogeneous image having Q = and = 1,1111. Then 

the parameters were estimated for the image by considering 
windows of size 4x4, 8x8, 16x16, 32x32 and 64x64. The mean 
Square error ^ ^ was also coirputed. This procedure was 
done for different window size and different ISEED values. 
From all such different error values we have taken the 
maximum error for a particular window size. As mentioned 
earlier the errors also vary randomly v/ith different ISEED 
values. The errors were then plotted against the windov/ 
size. As expected, the error was minimum for window size 
64x64, In Table 5.1, we present the error (maximum) for 
corresponding v^/indow size. 

Table 5.1 


VaNDOW 

ERROR (Max.) 

64 

0.9098x10”^ 

32 

0.2912x10”^ 

16 

0.3556x10“^ 

8 

0.9455x10“’^ 


4 


2 


0.293829 

0.352612 




X19 


The error vs. window size plot as shown in Fig. 5,6, 
suggests that there may be an exponential relation between 
error and window size. This experiment was carried out for 
homogeneous images with different parameter sets, I'ie have 
repeated the experiment by using different image sizes like 
48x48 etc. In all such cases, we have found that error and 
window size are related by the equation 

= Am"® (5.3) 

where , 


= mean square error 
w = window size 

For the image with 0 = [0.12 0.12]^ and -^ = 1.1111, if we 
estimate it as a 2 point problem then the eqn. (5.8) 
becomes 




1.456 


But estimating parameters for a 4 point neighbourhood for 
the same image gave 


£,,, = 1.46 


By changing either image size or the ISEED value it has 
been seen that A and B change. From simulation of different 
synthetic textures v^/e have seen that A and B depend on i) 
number of parameters, ii) image size and iii) the ISEED 
value. As the ISEED value can be varied randomly hence 
A and B also vary in a random fashion. For this reason 
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we have given different plots for this error, vs, window 
size in Fig. 5.7 to show the variation of error with ISHED 
and window size. Apparently, it is impossible to give any 
mathematical expression for A and B in terms of the above 
mentioned variables. We have computed the error of such 
a scheme for two cases 

i) Taking (9, 'p) as the feature vector 
ii) Taking 0 as the feature vector 


From the plots it is seen that the additional information 

obtained by including -Pin the feature vector does not change 

^ computed 

the error appreciably, VJe have also £ discrepency in 

terms of the w'indow size, where discrepency is defined as 
follows , 


If straight ■ 

line represents the actual segmen- 
ting line, then the hatched portion 

Fig. 5.8j Region of 

is defined to be the region of Confidence 

confidence as shown in Fig, 5.8 

Same concepts can be generalized to other segmenting curves. 
If the estimated segmenting line passes through this region 
we say that the segmentation procedure has converged 
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pn fvr: t..ly, ,in actual practice, the segmenting line 
m.sy a.-vi.sL.vi of this region. define discrepency 

as l.!u‘ raLlo ol number of blocks outside the region to the 
total numbor ol blocks. This discrepency follows (Fig. 5 . 9 ) 
an oxpon'-‘n tial relationship with the window size and it 
if rtaiur.t'tl by considering larger neighbourhood. One more 
i.bini) I. > here is that the discrepency-v^findow 

‘.j/u I'ulatiun iu independent of the segmenting curve. All 
ilif'uu uxpiT litHU! la. mentioned above were done for various 
l,ypi-f o! ran fmen tint! curves like st, line, circles, st, line 
with tloubln sIoih:- utc. and all the results for them 
wc? re con n i s t n n t . 

conclude our study of texture segmentation by 
dlG(.u'.sintj liotnc' of the disadvantages of the scheme 
df;vulopo<i here. 

For t!w purpose of elaboration, we consider the tv/o- 
rtrilon problem again. Let us assume that the two parameters 

fi) ^ 

ration. At the time of block estimation, we observe that 
estimate* ul a, or 0 . vary randomly about the actual or 

^ ^ in 


Q 


LC!l uc J 


>ay ti and ^ are the estimate of £1 and 02 * 


Li e ti ] 


and ; 


A 

"IL * ^ ^2H 

A 

^2L ^ 
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where 0^^ and are the lov/er value and upper value 

of the estimate of 9^ respectively. 

If we select 0^^, 0^ and the windov^ size such that 
the following conditions are satisfied^ 


i) I aiL - SsL I » 0 


ii) 

1 SlH 


1 » 

0 and 

iii) 

1 aiH 

- S2L i 

>> 

0 


where IX-Yj means the absolute difference between the vectors 
X and ^ , then the estimates of and ©^ are well separated 
and the clustering scheme will v\/ork perfectly. But suppose 
for a particular ©^ and ©^y the conditions i) and ii) are 
satisfied but instead of iii) we have 


a ^ o 
SlH 


then at the time of clustering, v/e may find intersections of 
such two clusters for some data points. Though we can 
assign these points to any of these clusters arbitrarily^ 
if the number of such points is large, then clustering 
scheme may give erroneous segmenting line. That is to 
say that if the differences between the corresponding 
elements in the vectors 9^ and ©2 is not greater than a 
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threshold value this scheme may fail to detect the actual 
segmenting line. This is because of the limitations of the 
estimation rule. If for some 0^^ and the scheme fails 
to detect the line properly we may say that statistically 
the texture is not a composite one and it may be treated as 
a homogeneous one. 





Actual 

FI6.5.5 THE ORIKNAl AND ESTIMATED VERSM OF 

different segmenting curve for a 
64x64 IMAGE. 
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FIG.5.6(o)ERROR VERSUS WINDOW PLOT FOR 
A 64 XS4 SYNTHETIC TEXTURE (4 
POINT ESTIMATION) 
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FI6.5.9 DISCREPENCY VERSUS WINDOW FOR 
A 64X64 synthetic IMAGE. 


chapter 6 

CONCLUSION 

In oux ^ "^^xtuxs modsllingji w© have considsrsci 

only the statistical raodeiiing concept and the special type 
of statistical considered here, is known as 'random 

field' model, in Random field models, we have considered 
both causal and nor'causai neighbourhoods. Procedures for 
implementation of ^uch models have been given in details. 

We have cited examples of application of such nodels to 
texture segmentaf^o’^ ^nd real image synthesis. A complete 
methodology has been given for frequency domain 

characterization ’’^^xtures. The application of 2D time 
series modelling "textural image has been studied for 
different Gauss ■*^®xtures. We have also studied the 
applicability an<^ ^^’^^letions of pseudo random arrays in 
image modelling# 

A new method fox image segmentation for a composite 
image by SAR is investigated. The area of image 

modelling is a vast in scope and the present study only 
considers one sp®^ific class of models. We have not specific 
class of models* have not even touched upon the entire 
field of syntac'b^*^ linage modelling. 
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Some future areas of investigation are cyclostationary 
models and time' series models other than linear Gaussian. 
Though linear non-Gauss ian time series models are very 
useful for modelling of some real data, the time series 
models, having random parameters are very useful for real 
life image modelling, A large amount of research into 
these models is now being concentrated on the construction 
and application of computationally efficient algorithms 
to determine order and obtain estimates of the unknown 
parameters which have desirable statistical properties [40] • 
The development of 2D counterpart of the varying parameter 
time series model may show a new direction for image 
modelling. 
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